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ABSTRACT 

In  this  thesis  we  derive  a  lattice  structure  to  realize  linear  phase  transfer 
functions  and  develop  an  adaptive  algorithm  for  continuously  updating  the  lattice 
reflection  coefficients.  The  lattice  structure  is  considered  because  of  its  superior  finite 
wordlength  performance  compared  to  transversal  structures.  The  adaptive  lattice 
algorithm  developed  in  this  thesis  has  been  applied  to  estimate  the  sinusoidal 
frequencies  as  part  of  Prony's  method.  Results  of  computer  simulation  supporting  the 
theory  are  reported. 
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I.  INTRODUCTION 

Every  finite-impulse  response  (FIR)  filter  xias  two  distinct  properties  [Ref.  5]  ; 
first,  it  is  always  stable;  second,  if  it  is  not  causal,  it  can  always  be  made  to  be  causal 
by  introducing  a  finite  delay.  FIR  filters  can  be  designed  so  that  their  frequency 
responses  have  an  exactly  linear  phase  characteristic.  FIR  filters  with  linear  phase  are 
important  in  applications  like  speech  processing  and  data  transmission,  because  in 
these  applications  a  nonlinear  phase  filter  is  harmful.  The  symmetry  property  of  the 
linear  phase  FIR  filters,  however,  helps  reduce  the  number  of  coefficients  by  nearly  one 
half  resulting  in  substantial  computational  savings. 

The  various  filter  realizations,  or  structures  that  are  frequently  considered  are  the 
direct  form,  the  cascade  form,  the  parallel  form,  and  the  lattice  form.  The  lattice  form 
realization  is  of  particular  interest  because  of  its  superior  numerical  performance,  and 
modularity  in  the  structure.  The  operation  of  a  lattice  filter  is  completely  described  by 
specifying  the  sequence  of  reflection  coefficients  that  characterize  the  individual  stages 
of  the  filter. 

Conventionally  an  adaptive  filter  is  composed  of  a  tapped  delay  line  or 
transversal  structure  with  adjustable  coefficients  or  weights  and  an  adaptive  algorithm 
which  updates  the  coefficients  continuously  based  on  some  performance  criterion.  The 
design  of  a  fixed  coefficient  filter  is  based  on  the  prior  knowledge  of  both  signal  and 
noise.  Adaptive  filters,  on  the  other  hand,  have  the  ability  to  adjust  their  own 
parameters  automatically,  and  their  design  requires  little,  or  no  a  priori  knowledge  of 
signal  or  noise  characteristics  [Ref.  14].  However,  the  designer  has  to  choose  the  order 
of  the  filter  and  the  type  of  the  algorithm.  Also  the  adaptive  filter  usually  requires  a 
large  initial  transient  time  (i.e.,  the  initial  filter  convergence  period). 

The  least-mean-square  (LMS)  adaptive  algorithm  minimizes  the  mean  square 
error  s(k)  by  recursively  altering  che  filter  coefficient  vector  3(k)  at  each  sampling 
instant.  The  original  Widrow-HofF  LMS  algorithm  is  B(k+i)  =  5(k)  +  2jle(k)X(k), 
where  fi  is  a  convergence  factor  controlling  stability  and  rate  of  adaptation  [Ref.  15]. 
The  algorithm  is  based  on  the  method  of  steepest  descent,  moving  b(k)  in  proportion 
to  the  instantaneous  gradient  estimate  of  the  mean  square  error.  The  successive 
orthogonalization  provided  by  the  lattice  offers  advantages  in  adaptive  convergence 
rate  which  cannot  be  achieved  with  tapped  delay  lines  [Ref.  17,18]. 
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Recently,  Prony's  method  has  been  applied  to  the  estimation  of  spectral  lines  in 
noise  [Ref.  11].  It  has  been  shown  that  the  Prony's  method  yields  better  spectral 
estimates  than  a  companion  spectral  estimation  technique,  called  Pisarenko's  method 
[Ref.  10].  Also  it  has  been  established  that  the  filter  structure  involved  in  Prony's 
method  has  a  iinear  phase  property  which  is  not  the  case  for  Pisarenko's  method  [Ref. 
11]. 

The  thesis  investigates  the  application  of  lattice  structures  in  Prony's  method  of 
spectral  line  estimation.  The  complete  solution  to  Prony's  method  consists  of  three 
steps:  (i)  representing  a  given  process  of  M  sinusoids  in  noise  in  terms  of  complex 
exponentials,  (ii)  finding  roots  of  a  symmetric  polynomial,  and  (iii)  estimating  the 
frequency,  phase  and  amplitude  information.  In  this  thesis,  however,  we  have 
emphasized  only  the  frequency  estimation  problem.  Here  we  derive  an  adaptive  lattice 
structure  to  realize  linear  phase  transfer  functions  which  will  be  used  to  estimate  the 
sinusoidal  frequencies  as  part  of  Prony's  method.  The  scope  of  the  thesis  consists  of 
obtaining  a  linear  phase  lattice  structure,  developing  a  least  mean  square  (LMS)  error 
based  adaptive  algorithm,  and  testing  the  lattice  structure  and  the  algorithm  by  means 
of  computer  simulation. 

The  thesis  is  organized  as  follows.  In  Chapter  II,  we  present  a  brief  review  of 
Prony's  method  of  representing  a  given  process  in  terms  of  a  set  of  complex 
exponentials,  and  then  address  the  problem  of  estimating  spectral  lines  using  this 
method.  We  show  that  the  original  Prony's  method  has  to  be  modified  slightly  in 
order  to  apply  it  to  the  spectral  line  estimation  problems.  In  Chapter  III,  we  discuss 
the  basic  concepts  of  the  linear  phase  FIR  filter  and  the  related  lattice  structure.  We 
show  three  examples  of  obtaining  the  lattice  structures  for  both  linear  phase  and  non- 
linear phase  FIR  transfer  functions  (Appendix  A).  In  Chapter  IV,  we  briefly  discuss 
the  least-mean-square  (LMS)  adaptive  algorithm  that  results  from  attempting  to 
minimize  the  mean  square  error  and  present  a  summary  of  some  LMS  algorithms  for 
lattice  that  have  been  reported  in  the  literature.  Also  included  are  some  computer 
simulation  results.  In  Chapter  V,  we  present  the  derivation  of  a  new  LMS  based  FIR 
adaptive  lattice  algorithm  and  extend  it  to  the  iinear  phase  case  as  well.  The  new 
algorithm  is  then  used  in  the  estimation  of  spectral  lines  in  white  noise.  Results  of 
simulation  are  included. 
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II.  PRONY'S  SPECTRAL  LINE  ESTIMATION 

A.  INTRODUCTION 

In  its  original  form,  Prony's  method  analyzes  processes  involving  simple,  or 
damped  harmonics  using  complex  exponential  functions  as  coordinate  functions.  On 
the  other  hand,  the  well  known  Fourier  analysis  consists  of  representing  a  given 
process  in  terms  of  a  set  of  sine  and  cosine  functions.  Recently,  Prony's  method  has 
been  applied  to  the  estimation  of  spectral  lines  in  noise  [Ref.  11].  It  has  been  shown 
that  the  Prony's  method  yields  better  spectral  estimates  than  a  companion  spectral 
estimation  technique,  called  Pisarenko's  method,  in  terms  of  bias,  spurious  responses, 
and  the  computational  complexity  [Ref.  10]. 

In  this  chapter,  we  present  a  brief  review  of  the  Prony's  method  of  representing  a 
given  process  in  terms  of  a  set  of  complex  exponentials,  and  then  address  the  problem 
of  estimating  spectral  lines  using  this  method.  We  show  that  original  Prony's  method 
has  to  be  modified  slightly  in  order  to  apply  it  to  the  spectral  line  estimation  problems. 

B.  PRONY'S  METHOD 

Consider  that  the  given  process,  x(k),  consists  of  n  distinct  sinusoids,  then  x(k) 
can  be  approximated  by  an  expression  of  the  form 


x(k)    =.  _L  [A  coscojk  +  Bj  sincOjk]  (2.1) 


where  the  co/s  are  sinusoidal  frequencies.   The  above  approximation  can  be  considered 
as  a  special  case  of  an  exponential  approximation  given  by  [Ref.  29] 


x(k)    =    L  C.  eJaik  (2.2) 

i=  1    ' 


where  m  =    2n  and  the  a.'s  are  identified  as  co.  and  -co..    The  values  of  co.  can  be 

i  it  i 

estimated,  by  Prony's  method,  assuming  that  the  data  are  known  at  k  =  0,  1,  ...,  N-l. 
From  Eq.(2.2),  we  can  obtain  the  following  set  of  equalities 
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c,    +  c2    +    ...     +  cm        =  x(0) 

C{  eJal     +  C2  eia2     4-    ...     4-  Cm  eJam     =  x(l) 

q  eJ2al         +  C2  eJ2a2       +    ...    +  Cm  ei2am         =  x(2)  (2.3) 


Ci  eJ(N.l)ai  +  ^  ej(N-Da2  +  ...  +  Cm  Q}(^-l)am  m  x(N.1} 

Now  we  have  a  set  of  N  equations  with  2m  unknowns,  namely,  C  and  a.  (i  =  1,  2,  ..., 
m)  which  can  be  solved  exactly  if  N  =  2m,  or  approximately  by  the  method  of  least 
squares  if  N>2m.  Also  note  that  the  N  equations  are  nonlinear  in  the  exponential 
terms  e^ai.    Let  e^ai,  i=  1,  2,  ...,  m,  be  the  roots  of  the  equation  [Ref.  29]. 


eJma  +  a1eKm-1)a  +  a2eKm-2)a  +  ...  +  a^ei*  +  am   =  0  (2.4) 

In  order  to  determine  the  coefficients  a.  (i  =  1,2,  ...,  m),  we  multiply  the  first  equation 

tVi 

in  Eq.(2.3)  by  O  ,  the  second  equation  by  «m.i,  •••,  the  m  equation  by  al  and  the 
(m+  1)  n  equation  by  1,  and  add  the  results.  In  this  way,  we  can  obtain  the  (N-m) 
linear  equations.   From  Eq.(2.2),  we  can  obtain  the  following  set  of  equalities 

x(m)  4-  x(m-l)a1  +  x(m-2)  a2  4-    ...    4-x(0)am    =    0 
x(m+l)  +  x(m)a:  4- x(m-l)  a2  4-    ...    +x(l)am   =    0 

(2.5) 


x(N-l)  4-  x(N-2)  a1  4-  x(N-3)  a2  4-    ...    +x(N-m-l)am    =    0 
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Since  the  data  x(k),  k  =  0,  1,  2,  ...,  N-l,  are  known,  this  set  of  equations  can  be  solved 
for  the  m  a's  if  N  2;  2m.  After  the  a's  are  determined,  e^ai,  i  =  1,  2,  ...,  m,  are  found 
from  Eq.(2.4).  The  set  of  equations  in  Eq.(2.3)  then  becomes  linear  in  terms  of  C,  i  = 
1,  2,  ...,  m.  The  C's  can  be  determined  from  the  first  m  equations,  or  determined 
approximately  by  applying  the  least  squares  method  to  the  entire  set  of  equations.  The 
nonlinearity  associated  in  finding  e^ai,  which  is  related  to  the  frequencies  jca  and  -jo).,  is 
concentrated  in  Eq.(2.4).   The  above  procedure  is  known  as  Prony's  method. 

Now,  in  the  case  of  Eq.(2.1),  since  the  roots  of  Eq.(2.4)  occur  in  reciprocal  pairs, 
Eq.(2.4)  must  be  invariant  under  the  substitution  of  e~Jai  for  e-lai,  so  that  we  must  have 
a2n=  l'  a(2n-l)  =  al'  -«  an+  1  =  an-l'   Thus  Ecl-(2-4)  becomes 


ej2nco  +  aj  ej(2n-l)(o  +  _  +  a^  ej(n+l)<fl  +  ^  ejnco 
+  on.!  J(n-1>»  +  ...  +  a{  ei«  +  i     =   o 

or 

ejnco  j  (ejnco  +  g-jna))  +  aj  (ej(n-l)co  +  e-j(n-l)a^  +  _  + 

anA  (ei03  +  e-i^)  +  an  ]     =     0 
since  e^nC0  #  0,  we  have 

2  cosnco  +  2  ot|  cos(n-l)o)  +■  ...  +  2  an_j  cosco  +  an   =    0  (2.6) 

Now  noting  that  m=2n,  and  applying  the  above  symmetry  of  a's  to  Eq.(2.5)  yields 
[Ref.  II] 
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(x(0)  +  x(2n)}  +  (x(l)  +  x(2n-l)}  a{  +  ...  +  {x(n-l)  +  x(n+l)}  on.j 
+  x(n)  an      =    0 

(x(l)  +  x(2n+l)}  +  {x(2)  +  x(2n)}  a{  +  ...  +  {x(n)  +  x(n+2)}  an-1 

+  x(n+l)an   =   0 

(2.7) 


(x(N-2n-l)  +  x(N-l)}  +  (x(N-2n)  +  x(N-2)}  a{  +  ...  +  {x(N-n-2) 
+  x(N-n)}  an_!  +  x(N-n-l)  an   =    0 


Eq.(2.7)  consists  of  a  set  of  N-2n  equations  and  n  unknowns.  This  set  can  be  solved 
exactly  for  the  a's  if  N=3n,  or  solved  approximately  by  the  least  squares  method  if 
N>  3n,  and  then  the  co's  are  determined  from  Eq.(2.6). 

C.       ESTIMATION  OF  FREQUENCY,  AMPLITUDE  AND  PHASE 

If  a  process  under  measurement  contains  an  unknown  number  of  sinusoids  of 
unknown  frequencies  and  amplitudes,  a  variant  of  Prony's  method  can  be  used  to 
determine  the  number  of  sinusoids  and  their  associated  frequencies  and  amplitudes.  As 
noted  above.  Prony's  method  is  a  technique  for  modeling  data  of  equally  spaced 
samples  by  i  linear  combination  of  exponentials.  An  exponential  curve  having  M 
exponentials  terms  can  be  determined  from  the  2M  data  measurements.  Each 
exponential  term  Aje  i  has  two  parameters  —  an  amplitude  A-  and  an  exponent  a; 
where  a.  can  be  real  or  imaginary.  For  the  case  where  only  an  approximate  fit  with  M 
exponentials  to  a  data  set  of  N  samples  is  desired,  such  that  N>2M,  a  least  squares 
estimation  procedure  is  used. 
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The  model  assumed  is  a  set  of  iM  exponentials  of  arbitrary  amplitude,  phase, 
frequency  and  damping  factor.  A  process  consisting  of  M  real  undamped  (a  is 
imaginary)  sinusoids  can  be  expressed  as 


x(k)  =  L1(azik  +  a.  Zj  k)  (2.8) 

1=1 

M 

=  L  A.  cos(2ixfkT  +  G.) 
i—  1  l  l 

with  Oj  =  (Ai/2)eJ  i  and  z.  =  e^27C  i   ,  where  A.  is  the  amplitude,  f{  is  the  frequency  and  Q.{ 

is  the  phase  of  the  i     sinusoid,  respectively,  and  T  is  the  sampling  interval. 

Finding  {A.,  0j,  fj  and  M  that  minimize  the  squared  error  is  a  difficult  nonlinear 

least  squares  problem.    Prony's  method  solves  two  sequential  sets  of  linear  equations 

with  an  intermediate  polynomial  rooting  step  that  concentrates  the  nonlinearity  of  the 

problem  in  the  polynomial  rooting  procedure  (similar  to  Eq.(2.4)). 

* 
Define  the  polynomial,  B(z),  which  has  z{  and  z{    as  its  roots,  given  by 

M  *         2M         <***  • 

B(z)  =  n  (z  -  z.)(z  -  z.  )  =    L    b.  z2M-J  =  0  (2.9) 

with  bQ=  1  and  the  b.  being  real  coefficients.   Since  the  roots  are  of  unit  modulus  occur 

e 
for  z.   Therefore,  Eq.(2.9)  can  be  written  as 


in  complex  conjugate  pairs,  then  Eq.(2.9)  must  be  invariant  under  the  substitution  z 


z2M  B(l/z)  =  z2M  LM  b.  zi-2M  =2Z    b.  zi  =  0  (2.10) 

j  =  0   J  j  =  0    J 

Comparing  Eqs.(2.9)  and  (2.10),  we  may  conclude  that  t>i =  ^2M-i  ior  i  =  ^'  •••  »^» 
with  Oq  =  b2^j  =  1.  Thus  the  requirement  for  complex  conjugate  root  pairs  of  unit 
modulus  is  implemented  by  constraining  the  polynomial  coefficients  to  be  symmetric 
about  the  center  element.  Hence  the  realization  of  B(z)  is  a  linear  phase  FIR  filter. 
Based  on  order  2M,  a  linear  prediction  error  can  be  written  as 
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£(k)  =  1     b  [x(k  +  j)  +  x(2M-k  +  j)]  (2.11) 

]  =  0  J 


which  reduces  the  number  of  coefficients  required  by  one-half. 

The  coefficients  b,.  b^,  ...   .  b      are  determined  in  a  least  sauares  fashion  bv 
i       _  m 

minimizing  the  total  squared  error. 


N-2M    o 
E  =    .|q      E2(k)  (2.12) 


which  yields  the  normal  equations 


m  N-2M 

=  o  k    i=0 


L   bk  [  L       [x(2M-k  +  i)  +  x(k  +  i)j  [x(2M-j  +  i)  +  x(j  +  i)]]  =  0  (2.13) 


for  k=  1 M 

This  equation  can  be  solved  recursively  for  increasing  order  M. 

From  the  estimated  {b.}  values,  the  {Zj}  are  determined  using  Eq.(2.9).   This  gives 
the  frequency  estimates. 


f  =  tan"1  [Im(zi)/Re(zi)]  /  2ttT  (2.14) 

To  obtain  the  {a}  a  second  set  of  normal  equations  is  solved, 


£     a.    t    zj  z.J   +   a.      1  z    1  z.  ] 
i=lL    i  -:  =  o  k-   l  v     i  =  ok 


=  .^ttRez^xQ)  (2.15) 

j  =  0  K 


fork=0,l,  ...  ,M 
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The  set  {a.}  provides  both  amplitude  (A.),or  power,  and  phase  (6.)  information: 

A.  =  |  a.  |  (2.16) 

6.  =  tan"1  [  Im(b.)/Re(b.)  ]  (2.17) 

In  the  foregoing  we  have  shown  that  the  frequency  estimation  as  part  of  Prony's 

method  requires  a  polynomial  with  a  linear  phase  property.  In  the  next  chapter,  we 
show  that  a  lattice  structure  can  be  utilized  to  implement  a  linear  phase  transfer 
function. 
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III.  LINEAR  PHASE  FIR  FILTERLNG 

A.  INTRODUCTION 

Every  finite-impulse  response  (FIR)  filter  has  two  distinct  properties  [Ref.  5]  ; 
first,  it  is  always  stable;  second,  if  it  is  not  causal,  it  can  always  be  made  to  be  causal 
by  introducing  a  finite  delay. 

FIR  filters  can  be  designed  so  that  their  frequency  responses  have  an  exactly 
linear  phase  characteristic.  FIR  filters  with  linear  phase  are  important  in  applications 
like  speech  processing  and  data  transmission,  because  in  these  applications  a  nonlinear 
phase  filter  is  harmful. 

The  impulse  response  sequence  of  a  linear  phase  FIR  filter  exhibits  a  kind  of 
symmetry,  e.g.,  h(n)  =  h(N-l-n)  for  n  =  0,  1,  ...,  (N/2)-l  (assuming  that  N  is  even). 
In  general,  FIR  filters  require  a  large  number  of  coefficients  to  adequately  meet  with 
the  required  filter  specifications.  The  symmetry  property  of  the  linear  phase  FIR 
filters,  however,  helps  reduce  the  number  of  coefficients  by  nearly  one  half  resulting  in 
substantial  computational  savings. 

The  various  filter  realizations,  or  structures  that  are  frequently  considered  are  the 
direct  form,  the  cascade  form,  the  parallel  form,  and  the  lattice  form.  The  lattice  form 
realization  is  of  particular  interest  because  of  its  superior  numerical  performance,  and 
modularity  in  the  structure.  Lattice  realizations  have  been  successfully  utilized  in 
filtering  and  spectral  analysis,  and  in  modeling  of  some  physical  process  like  speech, 
geophysical  data  etc.  [Ref.  8].  The  operation  of  a  lattice  filter  is  completely  described 
by  specifying  the  sequence  of  reflection  coefficients  that  characterize  the  individual 
stages  of  the  filter. 

In  this  chapter,  we  will  discuss  the  basic  concepts  of  the  linear  phase  FIR  filter 
and  the  related  lattice  structure.  And  finally,  we  will  show  three  examples  of  obtaining 
the  lattice  structures  realizing  for  both  phase  and  non-linear  nhase  FIR  transfer 
functions  (Appendix  A). 

B.  LINEAR  PHASE  FIR  FILTERS  [REF.  3] 

Let  (h(n)}  be  a  causal  finite  duration  sequence  defined  over  the  interval 
O^n^N-1  having  the  linear  phase  symmetry  property  given  by  h(n)  =  h(N-l-n).  The 
z  transform  of  (h(n)}  can  be  written  as 
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N-l 

H(z)  =  t  h(n)z"n  (3.1) 

n  =  0 


If  N  is  even,  then  £q.(3.1)  can  be  rewritten  as 

(N/2)-l  N-l 

H(z)=       ZJ     h(n)z"n+     t ...  h(n)z"n 

n  =  0  n  =  N/2 

Now  applying  the  symmetry  property  of  h(n)  in  the  second  term  on  the  right  hand  side 
yields 

H(z)  -W      h(n)z-«  ^iyi     h(N-l-n)z-(N-l-n) 

n=0  n=0 

which  can  be  simplified  to 


H(z)  =   (NL2>1  h(n)  {z"n  +  z<^-l-nh  (3.2) 

n  =  0 


If  N  is  odd,  then  Eq.(3.1)  becomes 


H(z)  _   KN-y/2M     h(n)[z-n  +  z-(N-l-n)]  +  h(M)z-[(N-l)/2]  (33) 

n  =  0  2 


Ifv/e  evaluate  Eqs.(3.2)  and  (3.3)  for  z=e^<0  ,  we  obtain  the  discrete  Fourier  transform 
of  the  filter  sequence  h(n)  defined  as 


H(ei°°)  =     *Y    h(n)  e-i0311  (3.4) 

n  =  0 


For  the  case  when  N  is  even,  we  then  have  the  discrete  Fourier  transform 
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H(e 


co)  =  e-jco((N-l)/2} 


(N/2)-l 
£       2h(n)cos[  o{n-(N-l)/2)}]1 

n  =  0 


(3.5) 


and  for  N  odd,  we  obtain 


H(ejW)-e-i»«N-l)/2][h((N.i)/2}  +  (N£)/2 

n  =  0 


cos[0){n-(N-l)/2}J].  (3.6) 


In  both  cases  above,  the  sums  in  brackets  are  real,  implying  a  linear  phase  shift 
corresponding  to  a  delay  of  (N-l)/2  samples.  Figures  3.1  and  3.2  show  the  direct  form 
realizations  of  an  FIR  filter  with  linear  phase  for  both  N  even  and  N  odd  cases. 


Figure  3.1     Direct-form  realization  for  an  FIR  filter  with  linear  phase  (N  =  even). 
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Figure  3.2     Direct-form  realization  for  an  FIR  filter  with  linear  phase  (N=  odd). 

C.       LATTICE  FILTERS  [REF.  4] 

The  basic  single  section  lattice  structure  is  shown  in  Figure  3.3  where  x(n)  is  the 
x(n)  and  fj(n)  and  gj(n)  are  generally  known  as  the  Forward  and  backward  prediction 
errors,  respectively.   The  defining  equations  of  the  lattice  are  given  by 


Wn> 

= 

Mn)  = 

x(n) 

tyn) 

= 

f0(n)  + 

^(n- 

1) 

gi(n) 

= 

kLf0(n) 

+  §o(n 

-1) 

(3.7) 


where  kj  is  called  the  lattice  reflection  coefficient.    If  another  section  is  cascaded  with 
the  first  one  the  resulting  equations  for  the  next  order  prediction  errors  are  given  by 
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f2(n)  =  fjCn)  +  k2gl(n-l) 
g2(n)  =  k/^n)  +  §1(n-l) 

Substituting  for  f,(n)  and  g^n-1)  from  Eq.(3.7)  yields 

f2(n)  =  x(n)  +  (kj  +  kjl^JxCn-l)  f  k2x(n-2)  (3.8) 

or 

f2(n)  =  x(n)  +  b2  ^(n-l)  +  b2  2x(n-2) 

and 

g2(n)=k2x(n)+{k1  +  k1k2}x(n-l)  +  x(n-2)  (3.9) 

or 

g2(n)  =  b2)2x(n)  +  b2  jx(n-l)  +  x(n-2) 

where  b2  j  =  kj(l  +  k2)  and  b2  2=k2< 

Thus,  by  induction,  for  a  cascade  connection  of  M  lattice  sections  we  have  the  general 
expressions 


M 

fM(n)=iJ:0    bM,ix(n"i) 

M 
gM(n)  =  JE     bM>M.iX(n-i)  (3.10) 


where  b^j  q=  1.  Eq.(3.10)  represents  the  FIR  filter  type  equations  to  obtain  the  Mm 
order  forward  and  backward  prediction  errors.  Now  taking  the  z-transform  of 
Eq.(3.10)  yields 
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M 


FM<Z>  -  is?0  *M,r  x(*) 


M 


GM(z)=i=L0  HlM-i2'1  X& 


(3.11) 
(3.12) 


For  an  unit  impulse  input,  i.e.,  when  X(z)=  1,  F^(z)  and  G^(z)  represent  the  forward 
and  backward  prediction  error  filter  transfer  functions,  respectively,  and 


GM(z)  =  z"MFM(z-1) 


(3.13) 


fm<k> 


Fm(z) 


G»(«) 


gm(K) 


.-i 


gm(k-l) 


m  +  i 

► 


m  + 


l(z) 


Wz' 


gm+i(^ 


Figure  3.3     Lattice  Section. 

By  inspection  of  Figure  3.3  we  see  that  the  M  stage  lattice  reflection  coefficient 
is  equal  to  the  FIR  filter  coefficient  b^  M  in  Eqs.(3.11)  and  (3.12).  In  other  words, 
we  have 
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kM=  bM,M  (3-14) 

In  fact,  the  FIR  like  prediction  error  filter  coefficients  can  be  iteratively  calculated 
starting  from  Eq.(3.14).  The  algorithm,  to  calculate  the  filter  coefficients  bvj  •,  also 
makes  use  of  the  well  known  lattice  property  that  an  M  order  lattice  contains  all 
prediction  filters  of  orders  m^M.  Now  consider  the  m  section  (1  ^m^M)  in  the 
cascade  connection  of  M  lattice  sections  which  can  be  described  by  the  following 
equations: 

Fm<z)  "  Fm-l<z>  +  kmz"lGm-l«  (3.15) 

Gm<z>  =  kmFm-l(z>  +    z''Gm.l(z>  (316) 

Eq.(3.16)  can  be  rewritten  as 

GmA(z)  =  zGm(z)  -  zkj^jCz)  (3.17) 

Substituting  Eq.(3.17)  into  Eq.(3.15)  yields 

Fm(z)=Fm-l(z)  +  WGm(z)-kmFm-l(z»  <3-18) 

Therefore,  the  (m-1)  n  order  forward  prediction  error  transfer  function  can  be  written 
in  terms  of  the  mu  order  forward  and  backward  prediction  error  transfer  functions  as 
follows: 

F-^fz)  -  k^G^fz) 

Fm-l(z)  -  J       g    m  (3-19) 

1   '  K  m 

where  km  =  bm     ¥*  1 .    Recalling  Eq.(3.13),  that  is, 

Gm(z)  =  z^F^z"1)  (3.20) 

By  substituting  Eq.(3.20)  into  Eq.(3.19)  we  find  that 
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F^fz)  -  k™z"mF™(  z"1  ) 
1   "   K  m 

Thus,  from  Eq.(3.21)  we  see  that  the  next  lower  order  polynomial  Fm_i(z)  can  be 
calculated  knowing  Fm(z).  Following  the  foregoing  procedure  we  can  find 
k.m_i  =  bm_j  m_j  from  Fm_i(z)>  and  also  obtain  Fm_2(z)-  This  way,  for  a  given  Mtjl 
order  polynomial  Fj^(z),  we  can  determine  all  the  lattice  reflection  coefficients  km, 
m=M,  M-l,  ...,  2,  1.  This  procedure  is  known  as  the  step-down  procedure  [Ref.  1], 
and  can  be  summarized  as  follows:   Let  the  given  M     order  polynomial  be 


FM(Z>  -  ;f0    bM,iZ"!  .  <3-22) 

and  by  replacing  M  with  m  we  have  an  mth  order  the  polynomial  for  m  lattice  sections 

FmW  =   jf0  bm,iz'1  <3'23> 

where  m=  M,  M-l,  ...,  2,  1.   We  now  define  another  polynomial  given  by 

As  we  step  through  the  procedure  from  m=  M  to  m=  1  Eq.(3.24)  can  be  expressed  as 

P   (Z-U  =   y1    b         -zm"^  (3  25) 

Define  km=  bmm  and  obtain  the  m-l"1  order  polynomial  as  follows: 
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Fm(z)  -  k™z"mF™  (z"1  ) 
Wz>  -       m\       g *^^  (3.26) 

1    "   K  m 


Substituting  Eqs.(3.23)  and  (3.24)  into  Eq.(3.26)  yields 


s 


>b    -z"1  -  k    7"mVb  -7m'1 

i=^Dmiz       Kmz  i:^QDm,m-iz 


i^)um-l,i 


k2 
K  m 


m  i  :_^rt  m,i  m:.^  m,m-i 

i'oVu2"1 <3-27> 

1      -      k2 

Then,  by  equating  the  coefficients  of  like  powers  of  z  on  both  sides  of  Eq.(3.27),  we 
obtain  the  computational  expression  for  the  (m-1)      order  polynomial  coefficients  as 


=     "mi "  kmbm.m-i 
1   -  k: 


v,  _       gu      m  m,m-i  /->  ->q\ 

bm-l,i ; — zr~^ ^3-28) 


m 

where  m  =   M,  M-1,  ...,  1  and  i  =   0,  1,  ...,  m-1,  k^  =  bmm  and  |km|  <  1  for  a 
minimum  phase  polynomial  Fm(z). 

D.       LINEAR  PHASE  FIR  LATTICE  FILTER 

In  the  foregoing  we  considered  obtaining  a  lattice  structure  from  a  given  FIR 
transfer  function,  and  vice  versa.  In  what  follows  we  will  deal  with  a  special  case  of 
FIR  filtering,  namely,  the  linear  phase  FIR  filter.  Let  {h(n)}  be  a  causal  finite  duration 
linear  phase  sequence  defmed  over  the  interval  O^n^N-1.    From  Eqs.(3.2)  and  (3.3), 

ihe  z  transform  of  (h(n)}  can  be  written  as 


HN.!(z)    -    FM(z)  +  z-D  GM(z)  (3.29) 
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Figure  3.4     Linear  Phase  FIR  Lattice  Filter. 

where  F-^(z)  and  Gyi(z)  are  forward  and  backward  filter  transfer  functions, 
respectively,  and  D  represents  the  number  of  unit  delays.  We  now  consider  the  N  odd 
and  N  even  cases  separately. 

iV  odd  case:   For  N  odd,  Eq.(3.29)  can  be  written  as 


HN_!(z)    -    a0  +  a^'1  +  a2z'2  +  ...  +  a(N_3)/2z 


•(N-3)/2    .    _  7-(N-l)/2 


+  a(N-3)/2z  +  - 


a->z 


.(N-3)  +  aiZ-(N-2)  ¥  „Z-(N-1) 


-1  -9 

aQ  +  a^z      +  a2z  *  +  ... 


,-(N-3)/2   .    7fi>o\n  7-(N-l)/2 

d(N-3)/2z  +  2  ^'2}  a(N-l)/2z 


+  a(N.3)/2z-(N+1)/2  +  ...  +  a2Z-(N-3)  +  ^-(N-2)  +  a()Z-(N-l) 
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Note  that  we  are   splitting  the  coefficient  a/^j^  mt0   two  coefficients   of  value 
0/2)a(N_1)/2  each. 

HN-l(z)    =    a0  +  ajz"1  +  a2z"2  +  ...  +  a(N_3)/2z"(N"3)/2  (3.30) 

+  (!/2)  a(N-l)/2z"(N"1)/2  +  z"(N"1)/2[(1/2)a(N-l)/2 
+  a(N-3)/2z_l  +  -  +  a2z"(N-5)/2  +  alZ-(N'3)/2 
+  a0z-(N-1)/2] 


HN-1(z)    =    FM(z)  +  z-W'W2  GM(z)  (3.31) 

where  the  order  of  the  polynomials  F\j(z)  and  Gvj(z),  M   =  (N-l)/2,  and  the 
number  of  unit  delays,  D  =  (N-l)/2,  the  forward  transfer  function, 


FM(z)    =    a0  +  a^1  +  ...  +  a(N.3)/2z-(N-3)/2  +  (i/2)eL(N_mz<^'^2  (3.32) 

FM(z-1)    =    a0  +  ajz  +  ...  +  a(N_3)/2z(N-3>/2  +  (l/2)a(N.1)/2z(N-1>/2  (3.33) 

and  the  backward  transfer  function, 

GM(z)    =    Z-(N-D/2  FM(z-';.  =  (l/2)a(N.1)/2  -  ^.3)/2zl  (3.34) 

+  ...  +  a,z-(N-3V2  +a0z-(N-'V2 
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N  even  case:   For  N  even,  Eq.(3.29)  can  be  written  as 

HN.i(z)    =    a0  +  a^1  +  a2z-2  +  ...  +  a(N.2)/2z"(N'2)/2  +  a(N-2)/2z"N/2 


+  ...  +  a2z"(N-3)  +  alZ-(N-2)  -  a0z-W) 


HN-1(z)    =    a0  +  ajz"1  +  a2z'2  +  ...  +  a(N.2)/2z-(N-2)/2  (3.35) 

+  ,-N/2  ra  4-       +  a  7-(N-6)/2 

+  z         la(N-2)/2  +  -  +  a2z 

+  alZ-(N-4)/2  +  a0z-(N"2)/2  ] 


HN-1(z)    =    FM(z)  +  z"N/2  GM(z)  (3.36) 

where  the  order  of  the  polynomials  F-yj(z)  and  G\Az),  M  =   {(N/2)-!},  and  the 
number  of  unit  delays,  D  =  N/2,  the  forward  transfer  function, 


FM(z)    -    a0  +  ajz"1  +  a2z"2  +  ...  +  a(N.2)/2z"(N-2)'2  (3.37) 

FM(z*1)    =    a0  +  aiz  +  a2z2  +  ...  +  a(N.2)/2z(N-2)/2  (3.38) 

and  the  backward  transfer  function. 

GM(z)    =    z-(N"2)/2  FM(z"1)  =  a(N.2)/2  +  ...  +  a2z-<N-6)''2  (3.39) 

+  ...  +  alZ-(N-4)/2+a0z-(N-2)/2 
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Figure  3.4  shows  a  linear  phase  FIR  lattice  realization  where  the  filter  output  is 

obtained  by  adding  the  M      order  forward  prediction  error,  and  the  M      order  delayed 

(by  D  unit  delays)  backward  prediction  error.    Combining  the  discussion  in  this  section 

and  that  in  Section  (III.C),  we  can  summarize  the  algorithm  for  converting  a  given 

FIR  transfer  function  into  a  lattice  structure  as  follows: 

(i)     Find  if  N  is  even,  or  odd 

(ii)   ForN  odd:  M  =  (N-l)/2,   D  =  (N-l)/2 

(iii)  For  N  even:  M  =  {(N/2)-l},   D  =  N/2 

(iv)  From  Fjyj(z),  obtain  M  reflection  coefficients  as  discussed  in  Eqs.(3.22)  to  (3.26) 

(v)   Implement  the  lattice  as  shown  in  Figure  3.4 

E.       LATTICE  REALIZATION  OF  A  GENERAL  FIR  TRANSFER  FUNCTION 

In  the  foregoing  we  have  considered  a  polynomial  Fjyj(z)  with  bjyj  q=  1.  However,  in 
general  we  have  b^Q^l.  In  this  section,  we  shall  modify  the  lattice  realization 
algorithm  presented  in  Section  (III.C)  to  suit  an  arbitrary  FIR  transfer  function  with 
DjV1  o  #  1  [Ref.  22].    The  m      order  polynomials    Fm(z)  and  Cfm(z)  can  be  obtained 

from  Eqs.(3.15)  and  (3.16): 

Fm«  =  smFm-l(z>  +  kmz"lGm-I<z>  <3-40> 

Gm(z)  =  kmFm-l(z)  +  sm  z"'Gm-l(z)  <3-41> 

where  the  coefficients  s„   =   b„  A  and  k„  =   b-,.  „  are  recognized  as  the  reflection 

m  m,u  m  m,m  ° 

coefficients.   Eq.(3.41)  can  be  rewritten  as 

Gm-l<z)  =  s'lm  z  Gm<z) "  s_1m  km  z  Fm-l<z)  ^3-42) 

Substituting  Eq.(3.42)  into  Eq.(3.40)  yields 

Fm<z)  =  *m  Fm-lW  +  ^m  WGmW  "  kmFm-l^  f3'43) 

Therefore,  the  (m-i)r  order  forward  prediction  error  transfer  function  can  oe 
written  in  terms  of  the  m  order  forward  and  backward  prediction  error  transfer 
functions  as  follows: 


31 


Fm-l<z>  = 


Figure  3.5     FIR  Lattice. 
sm  Fm<z) '  kmGm<z> 


b  m       K  m 


(3.44) 


where  sm#km.   Substituting   Eqs.(3.20).  (3.23)  and  (3.25)  into  Eq.(3.44)  yields 


m      m 


i^t) 


i^Vi/1  ~ 


s    f  h      .7"'  -  k    7"m?'b  -7m"i 


s2       -     k2 

s  m  ^  m 


UVl,iz     - 


m 


Y*  b„.  jz"1  -  k„  T*  b—  _  rz"1 
•^  q  m,i  m  -»-*_  q  m,m-i 


(3.45) 


s2        -      k2 
s  m  K  m 
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From    Eq.(3.45),   we    obtain   the    computational    expression    for    the   (m-l)th    order 
polynomial  coefficients  is 


u  sm  m.i  "    m^m.m-i  n  )A, 

bm-l,i II j3 "  ^-46> 

s  m  "  K  m 

where  m  -  M,  M-l,  ...,  1  and  i  =  0,  1,  ...,  m-1,  km  =  bm^m,  sm  =  bm0  and  sm>  1^ 

for  a  minimum  phase  polynomial  Fm(z).  It  may  be  noted  that  sm  =  1  for  m  =  1,  2, 
...,  M-l.  This  indicates  that  we  have  only  one  s-coefFicient,  i.e.,  Su,  which  requires  an 
extra  multiplication  in  the  M     lattice  section  as  shown  in  Figure  3.5. 
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IV.  LMS  ALGORITHM 

A.       INTRODUCTION      . 

Conventionally  an  adaptive  filter  is  composed  of  a  tapped  delay  line  or 
transversal  structure  with  adjustable  coefficients  or  weights  and  an  adaptive  algorithm 
which  updates  the  coefficients  continuously  based  on  some  performance  criterion. 

The  design  of  a  fixed  coefficient  filter  is  based  on  the  prior  knowledge  of  both 
signal  and  noise.  Adaptive  filters,  on  the  other  hand,  have  the  ability  to  adjust  their 
own  parameters  automatically,  and  their  design  requires  little,  or  no  a  priori  knowledge 
of  signal  or  noise  characteristics  [Ref.  14].  However,  the  designer  has  to  choose  the 
order  of  the  filter  and  the  type  of  the  algorithm.  Also  the  adaptive  filter  usually 
requires  a  large  initial  transient  time  (i.e.,  the  initial  filter  convergence  period). 

For  stationary  stochastic  inputs,  the  mean  square  error,  the  difference  between 
the  filter  output  and  an  externally  supplied  input  called  the  "desired  response",  is  a 
quadratic  function  of  the  filter  coefficients,  a  paraboloid  with  a  single  fixed  minimum 
point  that  can  be  sought  by  gradient  techniques  [Ref.  16]. 

In  the  previous  chapter  we  showed  that  the  operation  of  a  multistage  lattice  filter 
is  completely  described  by  specifying  the  sequence  of  reflection  coefficients  that 
characterize  the  individual  stages  of  the  filter.  In  this  chapter  we  briefly  discuss  the 
least-mean-square  (LMS)  adaptive  algorithm  that  results  from  attempting  to  minimize 
the  mean  square  error  and  present  a  summary  of  some  LiMS  algorithms  for  lattice  that 
have  been  reported  in  the  literature.  Also  included  are  the  computer  simulation  results. 

The  least-mean-square  (LMS)  adaptive  algorithm  minimizes  the  mean  square 
error  s(k)  by  recursively  altering  the  filter  coefficient  vector  B(k)  at  each  sampling 
instant.  The  original  Widrow-HofT  LMS  algorithm  is  B(k  +  l)  =  S(k)  +  2jie(k)X(k), 
where  \i  is  a  convergence  factor  controlling  stability  and  rate  of  adaptation  [Ref.  15]. 
The  algorithm  is  based  on  the  method  of  steepest  descent,  moving  b(k)  in  proportion 
to  the  instantaneous  gradient  estimate  of  the  mean  square  error. 

When  the  filter  input  is  stationary,  the  backward  prediction  errors  are  orthogonal 
to  each  other,  with  the  result  that  the  successive  stages  of  the  lattice  filter  are 
decoupled  from  each  other  [Ref  8].  This  means  that  the  global  optimization  of  a 
multistage  lattice  filter  may  indeed  be  accomplished  as  a  sequence  of  local  optimization 


problems,  one  at  each  stage  of  the  lattice  filter.  Accordingly,  it  is  a  straightforward 
matter  to  increase  the  order  of  the  lattice  filter  by  simply  adding  one  or  more  stages 
without  affecting  the  earlier  design  computations.  The  successive  orthogonalization 
provided  by  the  lattice  offers  advantages  in  adaptive  convergence  rate  which  cannot  be 
achieved  with  tapped  delay  lines  [Ref.  17,18]. 

B.       SUMMARY  OF  THE  LMS  ALGORITHM 

The  LMS  algorithm  uses  an  estimate  of  the  gradient  of  the  mean  square  error  obtained 
from  the  adaptive  linear  combiner  which  is  a  combination  of  a  transversal  structure 
and  an  adder.  The  adaptive  linear  combiner  can  be  shown  in  two  basic  ways, 
depending  on  whether  the  input  is  available  in  parallel  form  (multiple  inputs),  or  in 
serial  form  (single  input)  [Ref.  6]. 


d(k) 


y  (k) 


e  (k) 


xu-iiy-2(A  V^> 


Figure  4.1     Parallel  Form. 
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e  (k) 


Figure  4.2     Serial  Form. 

In  the  following  we  present  a  brief  derivation  of  the  LMS  algorithm.    Let  us  choose  the 
single  input  form,  then  the  filter  output  is  given  by 


y(k)  -  b0(k)x(k)  +  bjf^k-l)  +   ...   +  bN-1(k)x(k-N+l) 
=  £"  bi(k)  x(k-i) 

=  XT(k)  5(k) 


=  BT(k)X(k) 


(4.i; 


where  X  and  B  are  the  input  signal  vector  and  the  filter  coefficient  vector,  respectively 
and  (N-l)  is  the  order  of  the  filter.  The  input  signal  vector  X  and  the  filter  coefficient 
vector  B  are  defined  as 
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zoo  = 


x(k) 
x(k-l) 


x(k-N+l) 


B(k)    = 


bbf 


'N-lW 


/ 


Therefore,  the  output  y(k)  is  equal  to  the  inner  product  of  X(k)  and  B(k). 

The  error  e(k)  is  defined  as  the  difference  between  the  desired  response  d(k)  and 
the  actual  response  y(k), 


e(k)  =  d(k)  -  XT(k)S(k)  =  d(k)  -  ST(k)X(k) 


(4.2) 


The  purpose  of  the  adaptive  algorithm  is  to  adjust  the  filter  coefficients  of  the  adaptive 
linear  combiner  to  minimize  the  mean-square  error.  A  general  expression  for  mean 
square  error  as  a  function  of  the  filter  coefficient  values,  assuming  that  the  input 
signals  and  the  desired  response  are  statistically  stationary  and  that  the  filter 
coefficients  are  fixed,  can  be  derived  in  the  following  manner  [Ref  6].  The  squared 
error  is 


e2(k)  =  d2(k)  -  2d(k)XT(k)£(k)  +  BT(k)X(k)XT(k)B(k) 


(4.3) 


Taking  the  expected  value  of  both  sides  yields  the  mean  square  error, 


E[e2(k)]  =  E[d2(k)]  -  2  E[d(k)XT(k)]  g(k)  +  BT(k)  E[X(k)XT(k)]  B 


(4.4) 


Defining  the  vector  P  as  the  cross-correlation  between  the  desired  response  and  the 

input  vector,  we  have 


P   =    E[d(k)X(k)] 


(4.5) 


Similarly  the  input  autocorrelation  matrix  R  is  defined  as 
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R    =    E  [  X(k)  XT(k)  ] 


(4.6) 


Thus  the  mean-square  error  can  be  expressed  as 


£(k)  =  E[e2(k)j  =  E[d2(k)]  -  2  PT  g(k) .+  gT(k)  R  g(k) 


(4.7) 


The  gradient  Vs(k)  of  the  mean  square  error  function  is  obtained  by  differentiating 

Eq.(4.7)  with  respect  to  the  filter  coefficient  vector  as  follows: 


V£(k)  = 


=    -2  P  +  2  R  g(k) 


(4.8) 


The  LMS  algorithm  is  an  implementation  of  the  method  of  steepest  descent. 
According  to  this  method,  the  next  filter  coefficient  vector  is  equal  to  the  present  filter 
coefficient  vector  g(k)  plus  a  change  proportional  to  negative  of  the  gradient,  Ve(k): 


g(k+l)   =    g(k)-|i  7s(k) 


(4.9) 


The  parameter  \k  is  the  factor  that  controls  stability  and  rate  of  convergence.  In  other 
words,  the  first  term  on  the  right  hand  side  consists  of  the  past  information  and  the 
second  term  represents  the  new,  or  updated  information. 

The  LMS  algorithm  estimates  an  instantaneous  gradient  in  a  crude  but  efficient 
manner  by  assuming  that  e  (k),  the  square  of  a  single  error  sample,  is  an  estimate  of 
the  mean-square  error  and  by  differentiating  e^(k)  with  respect  to  3(k).  The  estimated 
gradient  is  given  by  the  following  expression 
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A 

Vs(k)  = 


3e2(k) 
5b0(k) 


de2{k) 


LfbN- 


(k) 


=    2  e(k) 


5e(k) 
5b0(k) 


5e(k) 


<3b 


(4.10) 


N-lW 


The  estimated  gradient  components  are  related  to  the  partial  derivatives  of  the 
instantaneous  error  with  respect  to  the  filter  coefficient  components.  Thus  the 
expression  for  the  gradient  estimate  can  be  simplified  to 


vs(k) 


■2  e(k)  X(k) 


(4.11) 


Using  this  estimate  in  place  of  the  true  gradient  in  Eq.(4.11)  yields  the  Widrow-Hoff 


LMS  algorithm 


B(k+1)    =    g(k)  +  2ne(k)X(k) 


(4.12) 


Since  the  filter  coefficient  changes  at  each  iteration  are  based  on  imperfect  gradient 
estimates,  we  would  expect  the  adaptive  process  to  be  noisy,  that  is,  it  would  not 
follow  the  true  line  of  steepest  descent  on  the  performance  surface.  The  LMS  algorithm 
can  be  implemented  in  a  practical  system  without  squaring,  averaging,  or 
differentiation  and  is  elegant  in  its  simplicity  and  efficiency.  Each  component  of  the 
gradient  vector  is  obtained  from  a  single  data  sample  without  perturbing  the  filter 
coefficient  vector. 

C.       ADAPTIVE  LATTICE  ALGORITHMS 

The  parameters  to  update  in  a  multistage  lattice  are  its  reflection  coefficients. 
Several  algorithms  have  been  proposed  in  the  literature  for  updating  the  lattice 
reflection  coefficients  [Ref.  7,8,17,18,20,23-28].  In  this  section  we  briefly  summarize 
some  of  those  algorithms.  Consider  a  lattice  filter  of  order  M.  For  stage  m  of  the 
filter,  the  flow  of  signals  is  described  by  the  following  pair  of  equations. 
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Figure  4.3    A  Single  Stage  of  Lattice  Filter. 


Uk)  ■■  fm.,(k)  +  kmf gm.,(k-i) 


nr 
gm(k) 


(4.13) 


where  m=  1,2,  ...  fM 

For  stage  m,  the  forward  reflection  coefficient  is  denoted  by  km  f  and  the  backward 
reflection  coefficient  is  denoted  by  km  ^.  fm(k)  and  gm(k)  are  the  forward  prediction 
error  and  backward  prediction  error  of  stage  m  respectively.    Before  presenting  the 

adaptive  algorithms  for  the  lattice,  let  us  quickly  summarize  some  non-adaptive 
reflection  coefficient  estimation  methods.  Later  on  we  can  obtain  the  adaptive 
updating  equations  as  approximations  to  these  methods. 
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1.  Non- adaptive  Methods 

When  the  lattice  coefficients  have  fixed,  non-adaptive  values,  several  methods 
have  been  proposed  for  computing  these  values  as  functions  of  the  correlation  statistics 
of  the  input.  Two  of  these  which  are  based  on  mean-square  error  (MSE)  minimization 
[Ref.28]  are  given  below. 

Method  1  :  In  this  method,  we  choose  kmf  to  minimize  the  mean  forward 
prediction  error  power,  E[f2m(k)]  and  1^  ^  to  minimize  the  mean  backward  prediction 
error  power,  E[g2m(k)].  Here  we  give  the  final  equations  for  km  r  and  km  j.  without 
going  into  the  details.   The  forward  and  backward  reflection  coefficients  are  given  by 


_    E[rm.i(fc)gm.,(k-i)] 

_     EiWOSm-l^-')] 


km,f '  p,„2 


where  both  km  c  and  km  -u  are  obtained  from  the  cross-correlations  between  the 
1)       order   forward 
prediction  error  powers. 


(m-l)in    order   forward    and    backward   prediction   errors    normalized   by    respective 


Method  2  :  For  a  single  channel  lattice  using  real  data  and  coefficients  we  can, 
however,  show  that  km£=km^  =  km.  Based  on  the  condition  that  km  f=  km  u  =  km, 
we  can  now  minimize  either  E[f2m(k)],  or  E[g2m(k)]  in  order  to  obtain  an  optimum  km< 
However,  it  seems  more  logical  to  minimize  the  sum  E[f  m(k)]  +  E[g  m(k)]  as  suggested 
in  [Ref.  28].   The  resulting  reflection  coefficient  equation  in  its  final  form  is  given  by 


„        _        2E[^n.1(k)gm,1(k-l)I 


Km.i 


E[f  m-lWl  +  Efg  m-l(k-l)] 


where  we  have  normalized  the  cross-correlation  term  with  respect  to  the  sum  of  mean 
prediction  error  powers.  In  both  of  these  methods,  the  coefficients  at  stage  m  can  be 
computed  independently  of  those  following  that  stage.    Thus,  optimum  values  can  be 
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computed  successively  along  the  structure  without  affecting  the  other  coefficients.   This 
phenomenon  lends  to  the  modular  nature  of  the  lattice  structure. 
2.  Adaptive  Methods 

An  instantaneous,  gradient-descent  adaptive  algorithm  minimizes  a  mean- 
square  error  criterion  can  be  derived  for  a  generalized  problem.  An  adaptive  lattice 
structure  is  suggested  by  incorporating  time  varying  coefficients  km^k)  and  km^(k)  in 
the  lattice  and  generating  algorithms  for  the  two  methods  described  above.  The 
resulting  procedures  are  : 

Method  1  :    Corresponding  to  method  1   of  non-adaptive  procedures  metioned 
above,  we  can  obtain  the  following  update  equations 


Wk+  !>  -  Wk>  +      g2  M        [«k)  gm-i(k-D] 

mA^  (4.16) 

kmfb»+  0  =  km,b«  +      g:    'U    (k)    fWk)  8mWl 

a  m-l,gw 

where  n  is  an  adaptive  step  size  parameter,  X  is  a  positive  weighting  constant,  and  the 
forward  power  estimate  at  the  (m-1)     stage,  <*2m_i  f(k),  is 

and  the  backward  power  estimate  at  the  (m-1)     stage,  G2m_\  u(k),  is 


Method  2  :    Eq.(4.15)  can  be  approximated  to  obtain  a  recursive  update  equation 

as  follows: 


k^k*  1)  =  km(k)  +  n  [f^g^k-l)  +  fmA(X)%m{k)]  (4.17) 
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For  the  basic  structure  of  the  single-channel  lattice,  ie,  km  e  =  km  ^  =  km,  the  reflection 
coefficients  or  filter  coefficients  km(k)  may  be  updated  using  a  modification  of  the 
LMS  algorithm  of  the  LMS  algorithm  [Ref.  17,20]. 


kjk+  1)  =  km(k)  +      /  [fm(k)gm.l(k-l)  +  f^OOg^k)]  (4.18) 

u  m-iw 

where  0~m_i(k)  is  the  power  estimate  at  the  (m-1)      stage.    Now  the  updated  power 
estimate  is 


<*2mA(k)  =  X  <*2m.ifc-l)  +  (i-X)  [%2mA(k-l)  +  f2mA(V]  (4.19) 

where  X  is  a  positive  weighting  constant  satisfying  the  criterion  0<X<  1  then  controls 
the  bandwidth  of  the  filter  and  the  resulting  power  averaging  time.  A  power  estimate 
is  required  at  each  stage  in  the  lattice  due  to  the  fact  that  the  forward  and  reverse  error 
sequences  have  decreased  power  with  increased  stage  number. 

Method  3  :  A  third  successful  method  of  implementing  an  adaptive  algorithm  for 
FIR  lattice  structures  has  been  reported  by  Griffiths  [Ref.  17,18].  This  algorithm  has 
been  originally  discussed  for  a  noise  canceller  application.  The  form  of  the  lattice 
noise-cancelling  adaptive  filter  is  shown  in  Figure  4.4.  The  lattice  noise-cancelling 
adaptive  filter  consists  of  an  M  stage  linear  prediction  lattice  for  the  reference  signal 
xr(k)  together  with  a  set  of  tap  coefficients  Vm(k)  which  provide  the  noise-cancelling 
subtraction  paths.  Griffiths'  algorithm  is  briefly  presented  in  the  following.  The 
update  equations  for  km(k)  are  given  by 


km(k  +  1)  =  km(k)  +         "  ffm(k)gm.1(k-l)  +  gm(k)fm-i(k)]  (4.20) 

where  the  power  estimate  at  the  (m-l)th  stage,  (72m_j(k)  is 


<r2m.i<k)  =  ^Vl^"1)  +  0-M  [g2m.i(k-l)  +  fViflO]  (4-21) 
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Figure  4.4    Lattice  Form  Implementation  of  Noise-cancelling  Filter. 

Each  coefficient  Vn(k)  can  be  determined  independently  of  Vn(k)  for  n>  m,  because  of 
orthogonalization  provided  by  the  lattice.   Thus  the  resulting  algorithm  is 


H 


Vm(k+1)=  Vm(k)  +  . 

7  mW 


W)  gni(k)J 


rm   '  &m 


(4.22) 


where  associated  power  measurement  y2    (k)  is 


m> 


Y2m(k)  =  XY2m(k-l)  +  (l-X)[gm(k-l)gm(k-l)] 


(4.23) 


and  em(k)  is  the  m     stage  error  signal  as  shown  in  Figure  4.4. 
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D.      SIMULATION  RESULTS 
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Figure  4.5    System  Identification  Experiment. 

For   the   purpose   of  computer    simulation,    let    us   consider   an   approximate 
algorithm  given  by 


km(k+l)  =  km(k)  +      -JL-    [gm,l(k-l)]e(k) 

a  mw 


(4.24) 


where 


<y2m(k)  =  X(T2m(k-l)  +  (1-X)  g2m.l(k-l) 


(4.25) 
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The  configuration  used  in  the  simulation  is  a  system  identification  experiment  as  shown 
in  Figure  4.5.   The  fixed  filter  transfer  function  considered  is  given  by 

H(z)  =  1  -0.89Z-1  +  0.25  z"2 
The  convergence  performance  of  the  LMS  algorithm  (as  given  by  Eq.(4.24))  can  be 
observed  by  plotting  the  error,  e(k).  versus  the  update  iteration,  k.  called  learning 
curves.  The  input  x(k)  to  both  fixed  and  adaptive  filters  is  a  white  noise  sequence  with 
zero  mean  and  unit  variance.  Figures  4.6  to  4.8  show  the  learning  curves  for  the  above 
example  where  we  have  used  three  different  values  for  the  adaptation  constant,  ft.  In 
the  next  chapter,  we  derive  a  new  algorithm  for  updating  the  reflection  coefficients, 
based  on  the  least  mean  square  principle  and  the  steepest  descent  algorithm. 
Improvement  in  convergence  speed  will  be  shown  using  the  new  algorithm. 
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Figure  4.6    Learning  Curve  (\i  =  0.01). 
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Figure  4.7     Learning  Curve  (ji  =  0.1). 


Figure  4.8     Learning  Curve  (jl  =  0.5). 
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V.  ADAPTIVE  LINEAR  PHASE  LATTICE  ALGORITHM 

A.  INTRODUCTION 

In  Chapter  III  we  dealt  with  the  realization  of  fixed  coefficient  FIR  lattice  filter  with 
both  linear  and  nonlinear  phase  characteristics.  We  reviewed  the  basics  of  LMS 
algorithm  and  some  adaptive  lattice  realizations  reported  in  the  literature  in  the 
previous  chapter.  The  three  adaptive  lattice  algorithms  discussed  are  direct 
approximations  of  their  non-adaptive  counterparts.  None  of  them  estimates  a  gradient 
as  required  by  the  LMS  algorithm.  In  this  chapter  we  present  a  derivation  for  a  new 
LMS  based  FIR  adaptive  lattice  algorithm  and  extend  it  to  the  linear  phase  case  as 
well.  The  new  algorithm  is  then  used  in  the  estimation  of  spectral  lines  in  white  noise. 
Results  of  simulation  are  included. 

B.  LMS  ALGORITHM  FOR  THE  FIR  LATTICE 

From  Eqs.(3.40)  and  (3.41)  the  FIR  lattice  equations  can  be  written  as 

yk)=smfm-l(k)  +  kmgm-l(k-1)  ^) 

Wk)  =  ^Sm-l^-1)  +  kmfm-l(k)  (5-2) 

where  m=  1,2,. ..,M,  sm=  1  for  m#M  and  km  are  the  lattice  reflection  coefficients.  A 
realization  of  Eqs.(5.1)  and  (5.2)  is  shown  in  Figure  5.1.  The  lattice  input  is 
x(k)  =  g0(k)  =  fQ(k).   The  output  of  the  filter  is 


y(k)  =  fM(k>  (5  3) 

=  s\l  fM-l^k)  +  kM  §M-l(k_l) 

Therefore,  the  error,  e(k),  is  given  by 


e(k)  =  d(k)  -  y(k) 

(5.4) 
=  d(k)-sMfM.1(k)-kMgM.1(k-l) 
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Figure  5.1     Adaptive  FIR  Lattice  Filter, 
where  d(k)  is  the  desired  signal.   The  objective  is  to  minimize  the  mean  square  error 


J  =  E  {  e2  (k)  } 


(5.5) 


The  gradient  of  J  with  respect  to  km  and  sm,  respectively,  are  given  by 


m    "     nv 


dJ 

lk~ 


p,. 


=  .1  Pf„A- 


E{e(k) 


Jv(k) 


5.6) 


and 
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=  -2  E{e(k)  fM.!(k)}  (5.7) 


^SM 


Then  the  LMS  algorithm  can  be  formulated  as  follows: 


kj(k+l)    =    kj(k)-|ik(-^-) 


3y(k) 

kj(k)  +  2  fik  E{e(k)  (  -|Li  )}  (5.8) 


J 


3y(k) 
=    kj(k)  +  2  Jik  e(k)  (  -^—  ) 


and 


sM(k+i)  -   sM(k)-^s(— -) 

=    sM(k)  +  2nsE{e(k)fM.1(k)}  (5.9) 

=    sM(k)  +  2  ns  e(k)  fM_j(k) 

where  ]i^  and  ns  are  the  adaptation  constant,  and  we  have  replaced  the  true  gradient 
by  its  instantaneous  estimate.   Defining 


z(k)   =    4~-  (5.10) 


yields 


kj(k+l)=  kj(k)  +  2^e(k)z(k) 
sM(k+  1)  =  sM(k)  +  2  fis  e(k)  fM-1(k) 
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where  j=  1,2,...,M.    The  next  step  is  to  estimate  the  gradient  vector  z(k).    For  this 
purpose,  define 


and 


Then  from  Eq.(5.3)  z(k)  is  given  by 


z(k)    =    <DMi(k) 

M'J  (5.14) 

-  ■mTOOm-IjW  +  W^M-lf1)  +  gM.l(k-])  6M,j 

where  §M  j  =  (3kj^/3k:).   Substituting  Eqs.(5.1)  and  (5.2)  to  Eqs.(5.12)  and  (5.13)  then 
we  have 


3f-  ,(k)  3g:  ,(k-l)  3k- 


3g:  i(k-l)  3f  i(k)  3k- 


therefore, 


<Difj(k)  =  sjCk^ijCk)  +  kjWT^jCk-l)  +  g^k-l)^  (5.15) 

Ti?j(k)  =  Si(k)¥Mj(k-i)  +  kjCk^^jCk)  +  q.iWy  (5.16) 


where 
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8i.i  ■  -%  (5-17) 

is  the  Kronecker  delta  function,  and  i=  1,2,  ...  ,M  and  j=l,2,  ...  ,M.    If  i=0,  then 

Eqs.(5.15)  and  (5.16)  are 


dx(k) 


4>oj(k)    =         ^  =    0  (5.18) 

,J  OK.: 


To,jW  =  °o,j(k)  =  °  <5-19) 

where  j=l,2,  ...  ,M.    Figure  5.2  shows  the  computation  of  gradient  elements  for  the 

adaptive  FIR  lattice  algorithm. 

From  the  Eqs.(5.15),  (5.16)  and  (5.19)  we  have 


*i,j(k>    "    Ti/k)    "    °  (5-20> 

where  l<i<j-l.   The  computations  of  gradient  elements  at  each  case  of  j=M,M-l,  ... 
,1  are  as  follows: 

Case  j-  M  : 

•o.m«  "  *o.mW  -  o 

<D;M{k)  =  SiikXD:.,  M(k)  +  k,(ki'F;.,  M(k-1>  +  gi.i(k-l)8iiM 

f  iiM(k)  =  s1(k)'Pi.1M(k-i)  +  w.,^)  +  q.!(k)8iiM 

where  i=  1,2,  ...  ,M.   From  Eq.(5.20),  every  gradient  element  equals  to  zero  except  final 
stage  elements  and  gradient  elements  of  the  final  stage  are 
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Vk-1)5Lj 


*ojuo 


^lUi 


¥0j(k) 


T  00 


Vkl5ij 


V2(k-1)5M-ij  Vi(Ic-1)5: 


^M-lj^ 


fM-2(k)5M-lj  fM-l°°    5M j 


Figure  5.2    Computation  of  Gradient  Elements  for  the  Adaptive  FIR  Lattice  Algorithm. 


Applying  the  Eqs.(5.17)  and  (5.20)  yields 
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*M,M«    "    fM-l«  (5-22) 

Case  j  =  M-l  : 

*0,M-l<k>    "    T0,M-1»    =    ° 

*i,M-l«  =  •iP^w.M-lft)  +  WW.lM  +  Si-lO'-^i.M-l 

1'i,M.,(k)  =  «i(k)THM.1(k.l)  +  kiW«'i.l,M-l«  +  fi-lW5i,M-l 

where  i=  1,2,  ...  ,M.   Applying  the  Eq.(5.20),  last  two  stage  elements  are  considerable, 
if  i=  M-l,  then  we  have 

*M-lfM-lW  "  »M-lW*M-2fM-lW  +  ^M-l^M-^M-l^1)  +  SM-^^M-LM-l 

TM-l,M-lW  "  sM-l(k)TM-2,M-l(k-1)  +  kM-l(k)°M-2,M-l(k)  +  fM-2^5M-l,M-l 
Applying  the  Eqs.(5.17)  and  (5.20)  yields 

®M-1,M-1«    =    gM-2(k-')  (5'23) 

*M-1.M-1«   =    fM-2<k)  <5-24) 

And  the  last  stage  terms  are 

^M.M-lW  =  sM(k)°M-l,M-l(k)  +  kMik^FM-l,M-lVk-i)  +  gM-l^'^M.M-l 

TM,M-l<k)  -  'M^M-l.M.l^1)  +  ^M^M-l.M-lW  +  fM-l(k)8M,M-l 
Using  the  Eqs.(5.17),  (5.23)  and  (5.24)  yields 
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*M,M-l(k>    =    sM(k)gM.2(k-l)  +  kM(k)  fM.2(k-l)  (5.25) 

TM,M-l(k)    =    sM(k)fM.2(k-l)  +  kM(k)gM.2(k-l)  (5.26) 

Finally, 
Case  j=  1  : 

*01(k)    -    T0)1(k)    =    0 


<Du(k)  =  s^k^^k)  +  k/k^^k-l)  +  gi.1(k-l)6i4  (5.27) 

Ti(1(k)  =  s^Ti.^^k-l)  +  ^(k^.^^k)  +  riA(k)ditl  (5.28) 

where  i=  1,2,  ...  ,M. 

The  LMS  algorithm  derived  in  the  foregoing  can  be  summarized  in  Table  1. 

Simulation  Results  : 

The  performance  of  the  algorithm  summarized  in  Table  1  has  been  observed  by 
computer  simulation.  The  configuration  used  is  the  system  identification  experiment  as 
discussed  in  Section  (IV.  D)  using  the  same  fixed  coefficient  filter.  The  learning  curves 
obtained  using  the  new  algorithm  are  shown  in  Figures  5.3  to  5.5.  Comparing  these 
with  the  learning  curves  in  Figures  4.6  to  4.8,  which  are  obtained  using  approximate 
algorithms,  we  observe  significantly  faster  convergence  rate  for  the  new  algorithm. 
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TABLE  1 
LMS  ALGORITHM  FOR  THE  FIR  LATTICE 

Initialization: 
K(0)    =  0 
sM(0)  =  1 

Lattice: 

x(k)  =  fQ(k)  =  gQ(k) 

fm(k)  =  ^(X)£Bl-1(k)  +  V^m-l^"1) 
gm(k)  =  sm(k)gm.l(k-l)  +  V^m-lW 
m  =  1,  2,    .  .  .  ,  M 

y(k)  =  %(k) 

Update  Equations: 

kj(k+l)  =  kj(k)  +  2  [fik/(T2k.(k)]  e(k)  Zj(k) 

sM(k+l)  =  sM(k)  +  2  [ys/(72s(k)]  e(k)  %.1(k) 

where  flj.  and  \is   are  the  adaptation  constants, 

(T2k.(k)  =  X  C72k.(k-1)  +  (1-X)  02M/j(k) 

and 

d2s(k)  =  X  (J2s(k-1)  +  (1-X)  f2M.1(k) 

are  estimations  of  power  in  z^(k)  and  fjy^-^k), 
respectively  and  X  is  a  positive  weighting  constant, 
0<X<1. 

Gradient    Vector  Elements: 

01#j(k)  =  si(k)0»i.lj(k)  +  ki(k)Ti.1/j(k-l) 

y±/j(k)  =  s-i(k)V1_1/j(k-l)  +  ki(k)<Di„1/j(k) 

+  fi-i(k)5i,j 
Zj(k)   =  0M/j(k) 
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Figure  5.3     Learning  Curve  (]i  =  0.1). 


Figure  5.4     Learning  Curve  (\i  =  0.3). 
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Figure  5.5     Learning  Curve  {]i  =  0.5). 


C.       LINEAR  PHASE  FIR  LATTICE  ALGORITHM 

We  now  extend  the  LMS  algorithm  derived  in  the  previous  section  to  the  linear 
phase  FIR  lattice  filter. 
From  Eq.(3.29),  output  of  the  linear  phase  FIR  lattice  can  be  written  as 


y(k)  =  fM(k)  +  gM(k-D) 


(5.29) 


The  lattice  input  is  x(k)  =sfQ( k)=gQ(k).    Substituting  Eqs.(5.1)  and  (5.2)  in  the  output 
equation  of  the  filter  yields 


y(k)  =  sM  flVM(k)  +  kM  gM.!(k-l)  +  sM  gM.1(k.D.1}  +  kM  fM.!(k-D|5 
"  SM  [fM-l<k)  +  8M-l(k-D-l)]  +  kM  lfM-l(k-D)  +  S.M-l^-1)! 

and  the  error,  e(k),  is  given  by 
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Figure  5.6    Adaptive  Linear  Phase  FIR  Lattice  Filter. 


e(k)  =  d(k)  -  y(k) 


"  d<k)  •  SM  trM-l<k)  +  gM-l(k-D-l)]  •  kM  lfM-l(k-D)  +  SM-1^-1) 


5.31) 


where  d(k)  is  the  desired  signal.  A  schematic  of  the  linear  phase  FIR  lattice  realization 
is  shown  in  Figure  5.6.  From  Eq.(5.31)  we  see  that  y(k)  is  a  function  of  both  s^j  and 
k^j.  But  frvj.j(k)  and  g^j(k-l)  are  functions  of  k^j_j  and  so  on.  Therefore,  in  order 
to  minimize  the  mean  square  error,  we  define  a  cost  function 


J  =  E  {  e-  (k)  } 


(5.32) 


and  a  gradient  element 
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5y(k) 

z(k)    =    -^-  (5.33) 

Then  following  the  treatment  in  Eqs.(5.6)-(5.11)  the  LMS  algorithm  for  the  linear 

phase  lattice  can  be  written  as 


kjCk+1)-  kj(k)  +  2nke(k)z(k) 

sM(k+l)  =  sM(k)  +  2  ms  e(k)  [fM.!(k)  +  gM-1(k-D-l)] 

where  j=  1,2,...,M.   From  Eqs.(5.29)  and  (5.30),  the  gradient  vector  z(k)  is  given  by 


z(k)   =   0M[j(k)  +  *MJ(M» 

+  TM.j/k-oi  +  !fVM(k-D)  -  gM.^k-ni  8Mj 

where  6^  ;  =  (dk^j/dkA  and  <Dm:  and  Tmj  are  obtained  on  the  same  lines  as  in 
Eqs.(5.15)-(5.17).   The  resulting  LMS  algorithm  is  summarized  in  Table  2. 

Simulation  Results  : 

The  performance  of  the  algorithm  summarized  in  Table  2  has  been  observed  by 
computer  simulation.  The  configuration  used  is  the  system  identification  experiment  as 
discussed  in  Section  (IV.  D).  We  have  used  two  different  fixed  coefficient  linear  phase 
FIR  filter  examples,  one  with  N  even  and  the  other  with  N  odd.   They  are 

H3(z)  =  0.154  +  0.462Z'1  +  0.462z'2  +  0.154z'3 
and 

H4(z)  =  0.15  -  0.45z_1  +  0.36z'2  -  0.45z"3  +  0.15z"4 
The  learning  curves  obtained  using  the  new  algorithm  are  shown  in  Figures  5.7  to  5.10. 
Figures   5.7   and   5.S    are   obtained  learning   curves   with   linear  phase    FIR   transfer 
function  H-j(z).    Figures  5.9  and  5.10  are  obtained  learning  curves  with  linear  phase 
FIR  transfer  function  H^(z). 
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TABLE  2 
LINEAR  PHASE  FIR  LATTICE  ALGORITHM 

Initialization: 

K(0)    =  0 

Sm(0)  =  1 
Lattice: 

x(k)  =  f0(k)  =  g0(k) 

fm(k)  =  sm(k)fm„l(k)  +  yklg^^k-l) 

gm(k)  =  sm(k)gm-1(k-l)  +  Vk)fm-l(k> 

m  =  1,  2,  .  .  .  ,  M 

y(k)  =  sM  [fM-1(k)  +  gM_1(k-D-l)]  +  kM  [fM-1(k-D) 

+  gM-1(k-l)] 
Update  Equations: 

k,(k+l)  =  k,(k)  +  2  [Jik/<r\  (k)]  e(k)  Zi(k) 

sM(k+l)  =  sM(k)  +  2  [fis/(J2s(k)]  e(k)  [f^^k) 

+  gM-i(k-D-i)] 

where  \i-^   and  ^s  are  the  adaptation  constants, 

(72k.(k)  =  X  <T2k.(k-l)  +  (1-X)  [a>M/j(k)  +¥M/j(k-D)] 

and 

(72s(k)  =  X  (J2s(k-1)  +  (1-X.)  [fM-!(k)  +  gM_1(k-D-l)] 

are  estimations  of  power  in  z^(k)  and 
[  %_]_(  k)+gM_1(  k-D-1)]  ,  respectively  and  X   is 
a  positive  weighting  constant,  0<X^1. 
Gradient   Vector  Elements: 


<E> 


0/j(k)   =  T0/j(k)   =  0 
<D,  ,(k)  =  B^k)^.,  ,(k)  +  iqCk)^.]  .(k-1) 

X  ,  j  —  X   -A.  ,  J  X  X—  /J 

+  gi_iCk-i)«i#j 

T^^k)  =  Si(k)Ti-1/;j(k-l)  +  k^kjO^^^k) 

+  fi-l(k)5i/j 
Zj(k)   =   *M/j<k>  +  TM,j^k"D) 
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Figure  5.7     Learning  Curve  (jn  -  0.1). 
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Figure  5.8     Learning  Curve  (ji  =  0.3). 
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Figure  5.9     Learning  Curve  (]i  -  0.03). 
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Figure  5.10     Learning  Curve  (ji  =  0.1). 
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D.      SPECTRAL  ESTIMATION 

In  this  section,  we  extend  the  linear  phase  FIR  adaptive  lattice  algorithm  derived 
in  the  previous  section  to  the  estimation  of  spectral  lines  in  white  noise.  The  spectral 
estimation  problem  is  some  what  different  from  the  system  identification  experiment 
that  we  have  been  considering  so  far.  In  a  typical  spectral  estimation  problem,  we  are 
given  a  time  series  and  we  need  to  estimate  its  spectrum.  A  suitable  configuration  that 
is  frequently  used  for  this  purpose  is  shown  in  Figure  5.11.  Comparing  this  with 
Figure  4.5,  we  notice  that  we  have  the  given  time  series  at  the  input  of  the  adaptive 
filter  and  the  filter  updates  its  coefficients  to  minimize  the  mean  square  output.  In 
doing  this,  we  assume  that  the  input  process  x(k)  has  been  generated  by  passing  a 
white  noise  sequence  through  a  filter,  say  I(z),  and  the  adaptive  filter  attempts  to 
realize  an  inverse  filter,  say  H(z)=  1/I(z).  Considering  that  the  adaptive  filter  has 
sufficient  degrees  of  freedom,  the  output  of  the  filter  will  be  a  white  noise  sequence. 

The  adaptive  algorithm  summarized  in  Table  2  can  be  used  in  spectral  estimation 
after  appropriate  modification.    In  the  present  case  we  minimize  the  cost  function, 

J  =  E{y2(k)} 

-■) 
rather  than  J  =  E{e^(k)}.   The  resulting  update  equations  can  be  shown  to  be 

kj(k+  1)  =  kj(k)  -  2  [Hk/<J2k.(k)]  e(k)  Zj(k)  (5.36) 

sM(k+  1)  =  sM(k)  -  2  Uis/c2s(k)]  e(k)  [fM4(k)  +  gM_1(k-D-l)]  (5.37) 

where  Z:(k),  <Jk.(k)  and  ffs(k)  are  as  defined  in  Table  2. 

Using  the  adaptive  algorithm  in  Eqs.(5.36)  and  (5.37)  we  can  obtain  values  of  the 
reflection    coefficients    and    the    s^    coefficient.     We    then    obtain    the    equivalent 

polynomial,  B(z),  by  going  through  the  standard  step  up  procedure  given  by  [Ref.  I]: 

bm,0    =    sm 

b  =    k 

m,m  m 

bm,i    =    smbm-l,i  +  kmbm-l,m-i 
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where   i=  l,2,...,m-l    and   m=  1,2,.. .,M.     Finally   the   required   linear   phase    transfer 
function  is  obtained  as  follows: 

HN-1<Z)    =    FM<Z)  +  z"°  GM<Z) 
where  Fjyj(z)=  B(z)  and  G-yj(z)=  z"MB(z_1).   The  spectrum  is  computed  as  follows: 


S(f)    = 


h0  +  h{  e-ilco  +  ...  +  hN_!  e-li*-1)**  |2 

where  (0  =  27tf.   and   f  is  the   normalized   frequency  (with   respect   to   the   sampling 
frequency)  in  the  range,  0  ^  f  <  0.5. 

Simulation  Results  : 

The  input  process  x(k)  consists  of  a  signal  in  noise,  given  by 
x(k)   =    s(k)  +  w(k) 
where  s(k)  may  be  a  single,  or  multiple  sinusoids  and  w(k)  is  a  zero  mean  unit  variance 
white  noise  sequence.    In  the  following  we  consider  several  examples  using  parameters 
ranging  from  a  single  sinusoid  to  4  sinusoids,  SNRs  from  30dB  to  lOdB,  and  filter 
order,  M,  from  2  to  30. 

Example  I:   We  consider  a  single  sinusoid  given  by 

x(k)    =    V2  cos(27tfk)  +  w(k)  (5.38) 

where  f=0.15,  SNR=30dB. 

Figures  5.12.  5.13,  5.14,  and  5.15  are  plots  of  frequency  spectrum  with  different 
order  of  lattice  Vi  and  adaptation  constant  \i. 

Figure  5.12  is  the  plot  of  the  case  M  =  2  and  ^  =  0.015.  We  see  that  the  peak  is 
at  f=0.15  and  no  spurious  responses. 

Figure  5.13  shows  the  frequency  spectrum  of  M  =  4  and  ^  =  0.01.  We  observe 
that  one  peak  is  at  f=  0.15  and  a  spurious  response  whose  magnitude  is  about  -53dB  at 
f=0.37. 
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Figure  5.11     Spectral  Estimation  Mode. 

Figure  5.14  shows  the  frequency  spectrum  of  M  =  10  and  ^i  =  0.03.  There  arc  one 
peak  at  f=0.15  and  four  spurious  responses  at  f=0.05,  0.25,  0.35  and  0.45.  The  largest 
spurious  response  is  about  -15  dB  at  f=0.45. 

Figure  5.15  shows  the  frequency  spectrum  of  M  =  20  and  j.1  =  0.03.  There  are  one 
peak  at  f=0.15  and  nine  spurious  responses.  The  largest  spurious  response  is  about 
-27  dB  at  f=  0.125. 

Through  the  Figures  5.12,  5.13,  5.14,  and  5.15  we  can  determine  that  the  best 
model  order  to  detect  the  one  sinusoidal  signal  is  VI  =  2. 

Next,  we  fix  the  order  of  lattice  M,  adaptation  constant  u  and  iteration  number 
k.  Figures  5.16,  5.17  and  5.18  are  the  plots  of  frequency  spectrum  with  changing  the 
signal  to  noise  ratio  (SNR). 

Figure  5.16  is  the  plot  of  M  =  20,  ^  =  0.05,  k=3000  and  SNR=30dB.  There  is  a 
peak  at  f=0.15  and  the  largest  spurious  response  is  about  -27dB  at  f=  0.125. 
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Figure  5.17  is  the  plot  of  M  =  20,  n  =  0.05,  k=3000  and  SNR=20dB.  The  largest 
spurious  response  at  f=  0.425  is  larger  than  the  desired  response.  At  this  case,  we 
cannot  estimate  the  frequency  of  sinusoid. 

Figure  5. IS  is  the  plot  of  M  =  20,  n  =  0.05,  k=3000  and  SNR  =  lOdB.  Three  of 
the  spurious  responses  are  larger  than  the  desired  response.  The  desired  response  is 
about  -13dB. 

From  Figures  5.16,  5.17  and  5.18,  if  the  SNR  is  getting  smaller  then  estimating 
sinusoidal  frequency  is  more  difficult. 

Example  2:   Next,  to  estimate  two  sinusoidal  frequencies  in  noise,  set  the  input  x(k)  as 

x(k)    =    V2  {  cos(27tf1k)  +  cos(27tf2k) }  +  w(k)  (5.39) 

where  the  normalized  frequencies  of  signals  fj  =  0.15  and  f2  =  0.25  and  set  SNR=30dB. 

Figure  5.19  shows  the  frequency  spectrum  of  M  =  4  and  ^1  =  0.02.  There  are  two 
peaks  at  fj  =  0.15  and  f2  =  0.25  and  no  spurious  responses. 

Figure  5.20  shows  the  frequency  spectrum  of  M  =  20  and  n  =  0.022.  There  are 
two  peaks  at  f|  =  0.15  and  f2  =  0.25,  and  eight  spurious  responses.  The  largest  spurious 
response  is  about  -20dB  at  f=  0.375. 

From  Figures  5.19  and  5.20,  the  best  model  order  to  estimate  the  frequencies  of 
two  sinusoidal  signals  is  M  =  4. 

Example  3:   Next,  to  estimate  four  sinusoidal  frequencies  in  noise,  set  the  input  x(k)  as 

x(k)    =    41.  {  cos(27tf1k)  +  cos(27if9k)  +  cos(27tf3k)  +  cos(27tf4k) }  +  w(k)  (5.40) 

where  the  normalized  frequencies  of  signals  fj  =  0.05,  f>  =  0.15,  f-  =  0.25,  '4  =  0.35.  and 
set  SNR=30dB. 

Figure  5.21  shows  the  frequency  spectrum  of  M  =  8  and  ]i  =  0.015.  There  are  four 
peaks  at  fj  =  0.07,  ^2  =  0.15,  fj  =  0.27  and  ^  =  0.375,  and  no  spurious  responses. 

Figure  5.22  shows  the  frequency  spectrum  of  M  =  30  and  ji  =  0.014.  There  are 
four  peaks  at  f^  =  0.05,  f2  =  0.15,  ^  =  0.25  and  f^  =  0.35,  and  eleven  spurious  responses. 
The  largest  spurious  response  is  about  -23dB  at  f=  0.45. 
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From  Figures  5.21  -and  5.22,  the  best  model  order  to  estimate  the  frequencies  of 
four  sinusoidal  signals  is  M  =  8. 
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Figure  5.12     Frequency  Spectrum  (1  sinusoid,  M  =  2,  ji  =  0.015). 
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Figure  5.13     Frequency  Spectrum  (1  sinusoid,  M  =  4,  ji  =  0.01). 
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Figure  5.14    Frequency  Spectrum  (I  sinusoid,  M  =  10,  n  =  0.03). 
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Figure  5.15     Frequency  Spectrum  (1  sinusoid,  M  =  20,  ^  =  0.03). 
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Figure  5.16     Frequency  Spectrum  (1  sinusoid,  M  =  20,  SNR=30dB). 
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Figure  5.17     Frequency  Spectrum  (1  sinusoid,  M  =  20,  SNR=20dB). 
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Figure  5. 1 8     Frequency  Spectrum  ( 1  sinusoid,  M  =  20,  SNR  =  lOdB). 


71 


o 

v—  - 

o 

1 

1 

MAGNITUDE(DB) 

-43  0       -29.0 

o 

i 

o 

0.00     0.05 

0.10 

3.15 

0.20         0.25         0.30 

FREQUENCY 

0.35 

0.40 

0.45 

Figure  5.19     Frequency  Spectrum  (2  sinusoids,  M  =  4,  ji  =  0.02). 
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Figure  5.20     Frequency  Spectrum  (2  sinusoids,  M  =  20,  ji  =  0.022). 
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Figure  5.21     Frequency  Spectrum  (4  sinusoids,  M  =  8,  \i  =  0.015). 
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Figure  5.22     Frequency  Spectrum  (4  sinusoids,  M  =  30,  \i  =  0.014). 
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E.       SUMMARY,  CONCLUSION  AND  SUGGESTED  FUTURE  WORK 

The  thesis  investigates  the  application  of  lattice  structures  in  Prony's  method  of 
spectral  line  estimation.  The  complete  solution  to  Prony's  method  consists  of  three 
steps:  (i)  representing  a  given  process  of  M  sinusoids  in  noise  in  terms  of  complex 
exponentials,  (ii)  finding  roots  oi  a  symmetric  polynomial,  and  (iii)  estimating  me 
frequency,  phase  and  amplitude  information.  In  this  thesis,  however,  we  have 
emphasized  only  the  frequency  estimation  problem. 

The  underlying  theory  of  Prony's  method  has  been  briefly  reviewed  in  Chapter  II. 
By  using  a  modified  Prony's  approach,  we  observed  that  the  frequency  estimation  (as 
part  of  Prony's  method)  requires  a  polynomial  with  a  linear  phase  property.  The 
basics  of  the  linear  phase  FIR  filter  and  the  lattice  structure  have  been  included  in 
Chapter  III.  Also  we  have  shown  that  a  linear  phase  FIR  transfer  function  can  be 
realized  from  an  FIR  lattice  using  D  number  of  additional  unit  delays  and  an  adder, 
where  D  is  determined  by  the  order  of  the  linear  phase  polynomial. 

The  principles  of  adaptive  filtering  and  the  LMS  algorithm  have  been  briefly 
studied  in  Chapter  IV.  We  have  derived  a  new  LMS  algorithm  for  the  FIR  non-linear 
phase  ans  linear  phase  transfer  functions  in  Chapter  V.  The  problem  spectral 
estimation  using  the  new  adaptive  algorithm  has  been  addressed,  and  the  simulation 
results  are  included. 

The  significant  outcome  of  the  work,  is  the  derivation  of  an  LMS  type  algorithm 
for  a  linear  phase  lattice  structure  and  its  application  to  spectral  line  estimation  as  part 
of  Prony's  method.  As  has  been  shown,  the  new  algorithm  yielded  faster  convergence 
rate  compared  to  previously  reported  approximate  algorithms.  The  application  of  this 
algorithm  to  spectral  estimation  produced  high  spectral  resolution  as  illustrated  by  the 
simulation  results. 

However,  we  have  observed  spurious  spectral  responses,  when  the  model  order  is 
higher  that  the  required.  Also  we  have  observed  that  the  SNR  performance  of  the 
algorithm  needs  to  be  improved.  For  input  SNRs  less  than  about  10  dB,  the  algorithm 
performed  poorly.  Finally,  we  have  dealt  with  only  the  frequency  estimation  part  of 
Prony's  method.  In  a  future  work,  some  of  the  above  mentioned  problems  may  be 
taken  up. 
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APPENDIX  A 

EXAMPLES  OF  OBTAINING  A  LATTICE  STRUCTURE  FOR  THE  FIR 

TRANSFER  FUNCTION 

Example  A.I 

Obtaining  a  lattice  structure  for  the  general  FIR  transfer  function.    The  unit 
sample  response  of  a  FIR  transfer  function  is  given  by 
H3(z)  =  0.5  +  0.25z_i  +  0.125z-2  +  0.0625z-3 

Solution:   Given  unit  sample  response  is 

H3(z)  =  0.5  +  0.25z_1  +  0.125z'2  +  0.0625z'3  (A.l) 

Using  Eq.(3.22),  we  have  polynomial  for  3  sections 

H3(z)-b3|0+    b3?1z-1+b3)2z-2    +b3)3z-3  (A.2) 

Comparing  Eqs.(A.l)  and  (A.2),  we  have 
b3)0  =  0.5 


b3,l  -  °-25 


b3>2  =  0.125 


(A.3) 


b3  3  =  0.0625 
Starting  with  m=  3  we  have  from  Eq.(3.46) 

k3  =  b3  3  =  0.0625 

s3  =  b30  =  0.5 

Now,  we  need  to  generate  the  coefficients  for  tWz)  and  from  Eq.(3.46) 
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(A.4) 


<;     h      •  -  k     h 

Vl,i 7~2 Z~T <A-5) 

And  for  m=  3  and  i  =  0  we  have 

h       _    s3b3.0  -  k3b3,3 
b2,0 ,2.k2 


•3-   -   k3 


(0.5X0.5)  -  (0.0625X0.0625) 
b?  0  =  7ZT7> 7TZ7T77) =   l  (A-6) 


2'°  "  (0.5)'2     -     (0.0625)2 


Next,  for  m  =  3  and  i  =  1  we  get 

u       _       s3b3.1  -  k3b3.2 
b2'l  ^~I^ 


(0.5X0.25) -(0.0625X0.125) 

tH  1  =  -j j =  0.47619  (A.7) 

2'*  (0.5r  -  (0.0625)2 

and  finally 

h       -       s3b3.2  "  k3b3.1 

b2,2 7*1 Z~t 

s3     "   k3 

(0.5X0.125) -(0.0625X0.25) 

b?  9  =  7 7 =  0.19048  (A. 8) 

l'L  (0.5)2  -    (0.0625)2 

from  Eq.(3.46),  we  have 

k-2  =  b2<2  =  0.19048 


s2  =  b2>0  =  1 


The  new  polynomial  is 

H2(z)  =  b2)0  +  b2Jz-1  +  b2)2z- 


2 


(A.9) 
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H2(z)  =  1  +  0.47619z_1  +  0.19048Z'2  (A.10) 


for  m=  2  and  i  =  0,  we  have 

h  b2.0  '  k2b2.2 

bU  1    -   I? 


1      -(0.19048)(0.19048) 

=  1  (A.ll) 


"1,1  =  - 

i 

-    (0.19048)2 

and  finally 

h*   ,   —  - 

b2,l 

-  k^b-)  | 

bU 

1    - 

V 

0.47619 -(0.19048)(0.47619) 

b,  ,   = -n i_  =  0.40000  (A.12) 

l'x  1     -     (0.19048)- 


From  Eq.(3.46),  we  have 
kl  =  bl  1  =  °-40000 


(A.13) 


sl  =  bl,0  =  l 


Therefore,  reflection  coefficients  and  s  coefficients  of  the  FIR  lattice  structure  are 
kj  =    0.40000 


k2 

= 

0.19048 

k3 

= 

0.0625 

sl 

= 

1 

S-7 

= 

1 

s3 

= 

0.5 
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Example  A. 2  (N=  4:even) 

Obtaining  a  lattice  structure  for  the  linear  phase  FIR  transfer  function.   A  linear 
phase  FIR  transfer  function  is 

H3(z)  =  0.154  4-  0.462z_1  +  0.462z"2  4-  0.1 54z"3  (A.  14) 
Solution:   Eq.(A.14)  can  be  written  in  the  form  of  Eq.(3.35) 

H3U)  =  3.q  +  ajz"1  +  z""  (a|  4-  a^z"")  (A.  15) 
Comparing  Eqs.(A.14)  and  (A.  15),  we  get  the  following  relationships. 


a0  -  0.154 


a  {  =  0.462 


(A.16) 


To  determine  the  lattice  reflection  coefficients  of  the  corresponding  linear  phase  FIR 
lattice  structure  and  the  number  of  unit  delay,  we  use  Eq.(3.29)  which  is 

HN_!(z)    =    FM(z)  +  zD  GM(z)  (A.17) 

From  Eq.(A.14),  we  have  N  =  4(even),  the  order  of  the  polynomials  F^z)  and  Gjyj(z), 
M  =  (N/2)-l  =  (4/2)- 1  =  1,  and  the  number  of  unit  delays,  D  =  N/2  =  2.  Now, 
from  Eq.(3.37),  we  may  have  the  forward  prediction  error  transfer  function  as  follows: 


Fj(z)    =    a0  +  aj  z"1 


=    0.154    +  0.462  z_1 


for  .VI  =  1,  Eq.(A.iS)  can  be  rewritten  as 


(A.18) 


Fi<z)    =    b10  4-  bu  z"1  (A.19) 
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From  Eqs.(A.18)  and  (A.  19),  we  get  the  desired  reflection  coefficient  and  s  coefficient 
of  the  linear  phase  FIR  lattice  structure  as  follows: 

kl    =    bl,l    =    0A62 
sl    =   bl,0   =   °-154 

Example  A. 3  (N=  5:odd) 

Obtaining  a  lattice  structure  for  the  linear  ^phase  FIR  transfer  function.   A  linear 
phase  FIR  transfer  function  is 

H4(z)  =  0.15  -  0.45z_1  +  0.36z-2  -  0.45z-3  +  0.15z"4  (A.20) 

Solution:   Eq.(A.20)  can  be  written  in  the  form  of  Eq.(3.30) 

H4(z)  =  a0  +  ajz'1  +  (l/2)a2z-2  4-  z"2  {(l/2)a2  +  ajz"1  +  a0z*2}  (A.21) 

Comparing  Eqs.(A.20)  and  (A.21),  we  get  the  following  relationships. 

aQ  =  0.15 

aj  =  -0.45  (A.22) 

a2  =  0.36 

To  determine  the  lattice  reflection  coefficients  of  the  corresponding  linear  phase  FIR 
lattice  structure  and  the  number  of  unit  delay,  we  use  Eq.(3.29)  which  is 

HN-1(z)    =    FM(z)  +  zD  GM(z)  (A.23) 

From  Eq.(A.20),  we  have  N  =  5(odd),  the  order  of  the  polynomials  Fpvj(z)  and  G^(z), 
M  =  (N-l)/2  =  (5-1V2  =  2.  and  the  number  of  unit  delays,  D  =  (N-1V2  =  2.  Now, 
from  Eq.(3.32),  we  may  have  the  forward  prediction  error  transfer  function  as  follows: 

F2(z)    -    a0  +  Hzl  +  (1/2)  a2z"2 

(A.24) 
-    0.15   -0.45Z'1  +  0.18  z"2 

for  M  =  2,  Eq.(A.24)  can  be  rewritten  as 
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F2(z)    =    b20  +  b2>1  z'1  +  b2)2z'2  (A.25) 

From  Eqs.(A.24)  and  (A.25),  we  have 


k2    =    b22    =    0.18 

s2    =    b2)0    =    0.15 

Now,  we  need  to  generate  the  coefficients  for  Fj(z)  and  from  Eq.(3.46) 


s2b20-k?b2?         (0.15X0.15) -(0.18X0.18) 
\0  =         "  2'    .  r 2  '"     -      0.152-0.182 =  !  (A,26) 


.      ^2  1-^2,1      =(0.15X-0.45)-(0.18X-0.45) 
1,1  s22    -  k22  0.152-  0.182 

From  Eqs.(A.26)  and  (A. 27),  we  have 

kl    =    bl  1    =  -1-36364 

sl    "    bl,0    =    * 
Therefore,  reflection  coefficients  and  s  coefficients  of  the  linear  phase  FIR  lattice 
structure  are 


h 

= 

-1.363 
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h 

= 

0.18 

si 

= 

1 

s2 

= 

0.15 
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APPENDIX  B 
COMPUTER  PROGRAMS 

*  This  program  is  designed  for  the  system  identification  experiment.    * 

*  which  is  shown  in  Section  (IV. D).  The  learning  curves  can  be  obtained* 

*  by  plotting  the  error,  E,  versus  the  update  iteration,  k.  * 
************************************** 

* 

Variable  Definition 


ISEED 

AMU 

k 

M 

NB 

A(I) 

B(I) 

F(I) 

G(T) 

GD(I) 

SGL(I) 

W(K) 

YF(K) 

YL(K) 

ER(K) 


seed  for  the  random  number  generation  (white  noise) 

adaptation  constant 

time  index 

order  of  the  FIR  transfer  function  or  total  number  of 

lattice  sections 

number  of  iterations 

coefficients  of  the  FIR  transfer  function 

reflection  coefficients  of  the  lattice 

forward  prediction  error 

backward  prediction  error 

delayed  backward  prediction  error 

estimations  of  power 

input  of  both  FIR  and  lattice 

output  of  the  FIR  filter 

output  of  the  lattice  filter 

squared  error 

Variable  Declaration 


INTEGER  ISEED,K,I,J,N,M,NB 

REAL  A(100) ,B(100) ,F(100) ,G(100) ,GD(100) ,YF(10000) ,YL(10000) 

REAL  X(100) ,SGL(100) ,ER(10000) ,W(10000) ,Y(10000) ,AMU,E 

Set  Adaptation  Constant  |X  and  Number  of  Iterations  MB 


WRITE (5,1) 

FORMAT ( 5X , ' AMU ' , 5X , ' NB ■ ) 

READ(5,*)  AMU,NB 

Initialization 

DO  10  K=l,100 
A(K)=0 
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B(K)=0 

X(K)=0 

F(K)=0 

G(K)=0 

GD(K)=0 

SGL(k)=1.0 
10   continue 

DO  15  K=l, 10000 

YF(K)=0 

YL(K)=0 

ER(K)=0 

W(K)=0 

Y(k)=0 
15    CONTINUE 
E=0. 
R=l 

ISEED=343169 
M=2 

A(l)  ~      1.0 
A(2)  =  -0.89 
A(3)  =  +0.25 
* 

*  Random  Number  Generation 

(mean  zero,  unit  variance,  white  sequence) 
* 

DO   20   N=1,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 

W(N)  =  RN 

20    CONTINUE 
* 

*  FIR  Filter  Calculation 
* 

DO   30  K  =  1,NB 
X(1)=W(K) 

YF(K)=  A(1)*X(1) 
DO  40  1=1, M 

YF(K)=YF(K)+A(I+1)*X(I+1) 
40       CONTINUE 

DO  45  1=1, M 

X(M+2-I)=X(M+l-I) 
45       CONTINUE 
30   CONTINUE 

Lattice  Filter  Calculation 
* 

DO   50  K  =1,NB 
F(1)=W(K) 
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G(1)=W(K) 

DO   60   I  ■  1,M 

F(I+1)=F(I)+B(I)*GD(I) 

G(I+1)=GD(I)+B(I)*F(I) 

60       CONTINUE 

YL(K)=F(M+1) 
* 

*  Calculating  the  Error 

* 


E  =  YF(K)  -  YL(K) 


* 


*  Updating  the  Reflection  Coefficients 
* 

CALL  LMS  (B,GD,E,AMU,SGL,M) 
* 

ER(K)=E**2 
DO  70  J=1,M 
GD(J)=G(J) 
70       CONTINUE 

IF  (K.EQ.R)  THEN 

WRITE(6,300)  K,  B(l) ,3(2) ,ER(K) 
R=R+50 
END  IF 
Y(k)=E 
50   CONTINUE 
300   FORMAT(3X,I6,5X,3(F10.7,4X)) 

*  Plotting  the  Learning  Curve 
* 

CALL  PLOT(Y,N) 

STOP 

END 

* 

SUBROUTINE  LMS(B,GD,E ,AMU, SGL,N) 

REAL  B ( 100 ) , GD ( 100 ) , SGL ( 100 ) , E , AMU 
INTEGER  M 

DO  200  1=1 ,M 

SGL(I)=0.9*SGL(I)+0.1*(GD(I)*GD(I)) 
200      CONTINUE 

DO  210  1=1, M 

B(I)=B(I)+(AMU/SGL(I))*E*GD(I) 
210      CONTINUE 
RETURN 
END 
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SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N),X(10000) 

ISTP=N/10 

DO  10  J=1,N 
10   X(J)=J 
c    CALL  TEK618 
C    CALL  PRTPLT(72,6) 

CALL  SHERPA( 'PARKPARK' , 'A' ,3) 

CALL  ?HYSOR(l. ,1.) 

CALL  RE SET ('ALL') 

CALL  PAGE(8. 5,11.0) 

CALL  HWROT('AUTO') 

CALL  XINTAX 

CALL  AREA2D(5.0,2.8) 

CALL  HEIGHT (0.10) 

CALL  COMPLX 

CALL  SHDCHR(90. 0,1 ,0.002,1) 

CALL  HEADINC LEARNING  CURVE$ ', 100,2.0 ,1) 

CALL  XNAME ( ' ITERATIONS S ' , 100 ) 

CALL  YNAME ( ' ERROR$ ' ,100) 

CALL  MESSAGC ADAPTIVE  LATTICE(AMU=. 5 ,FIG4.8)$ ' , 100 ,3 .0, -0 .8) 

CALL  FRAME 

CALL  GRAF(0,ISTP,N,-3.0,1.5,3.0) 
c    CALL  THKCRV(0.02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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******pr0gram  2********************************************************** 

*  This  program  is  designed  for  the  system  identification  experiment     * 

*  using  the  LMS  algorithm  which  was  derived  in  Section  (V.B) .  * 

* 


* 


INTEGER  ISEED,K,I,J,N-,M,NB 

REAL  A( 100 ),B( 100) ,F( 100 ),G( 100) ,GD(100) ,x(100) ,YF( 5000), YL(5000) 

REAL  ER(5000) ,W(5000) ,SGL(100) ,PH(100/100) ,PS(100,100) 

REAL  PSD(100,100) ,GR(100) ,Y(5000) ,AMU,E 


*  Set  Adaptation  Constant  H   and  Number  of  Iterations  NB 

WRITE (5,1) 
1    FORMAT(5X, 'AMU' ,5X, 'NB') 
READ(5,*)  AMU,NB 

*  Initialization 
* 

E=0. 
R=l 

DO  5  K=  1,5000 
W(K)=0 
YF(K)=0 
YL(K)=0 
ER(K)=0 
Y(K)=0 
5    CONTINUE 

DO  10  K=l,100 
A(K)=0 
B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
SGL(K)=1.0 
GR(K)=0 
10       CONTINUE 

DO  15  K=l,100 
DO  15  L=l,100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
15       CONTINUE 
ISEED=343169 
M=2 
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A(l)  =   1.0 

A(2)  =  -0.89 

A(3)  =  +0.25 
* 

*  Random  Number  Generation 

*  (mean  zero,  unit  variance,  white  sequence) 

DO   20  -N=1,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 
W(N)  =  RN 

20       CONTINUE 

*  FIR  Filter 

DO   30  K  =  1,N3 
X(1)=W(K) 
YF(K)=  A(1)*X(1) 
DO  40  1=1, M 

YF(K)=YF(K)+A(I+1)*X(I+1) 
40  CONTINUE 

DO  45  1=1, M 

X(M+2-I)=X(M+l-I) 
45  CONTINUE 

30      CONTINUE 

*  Lattice  Filter 
* 

DO  50  K  =1,NB 
F(1)=W(K) 
G(1)=W(K) 
DO   60   I  =  1,M 

F(I+1)=F(I)+B(I)*GD(I) 
G(I+1)=GD(I)+B(I)*F(I) 
60  CONTINUE 

YL(K)=F(M+1) 
E  =  YF(X)  -  YL(K) 
* 

UDdatina  the  Reflection  Coefficients 


X 


CALL  MLMS  (PH,PS,PSD,GR,F,G,GD,B,SGL,E,AMU,M) 
ER(K)=E**2 
DO  70  J=1,M 
GD(J)=G(J) 
70  CONTINUE 

IF  (K.EQ.R)  THEN 

WRITE(6,300)  K,  B(l) ,B(2) ,B(3) ,B(4) ,B(5) ,B(6) ,B(7) ,B(8) ,E 

R=R+30 
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END  IF 

Y(K)=E 

c    WRITE(6,100)  K,W(K),YF(K),YL(K),E,ER(K) 

50      CONTINUE 

300   FORMAT(1X,I6,2X,9(F10.7,2X)) 

100   FORMAT(3X,I3,2X,5(F10.7,2X)) 
* 

*  Plotting  the  Learning  Curve 

* 

CALL  PL0T(Y,N) 

STOP 

END 
* 

* 

SUBROUTINE  MLMS (PH,PS ,PSD ,GR,F ,B ,BD ,R,SGL,EK, AMU,N) 

DIMENSION  PH(100,100) ,PS(100,100) ,PSD(100 , 100) ,F(100) ,B(100) 
1,BD(100),R(100),GR(100),SGL(101) 

GR(N)=BD(N) 

DO  200  1=2, N 

PH(I,1)=BD(N+1-I) 

PS(I,1)=F(N+1-I) 

IT=I-1 

DO  10  K=1,IT 

PH(I,K+1)=PH(I,K)+R(N+1-I+K)*PSD(I,K) 
10   PS(I,K+1)=PSD(I,K)+R(N+1-I+K)*PH(I,K) 

DO  20  K=1,IT 
20   PSD(I,K)=PS(I,K) 
200   CONTINUE 

DO  210  K=2,N 

210  GR(N+1-K)=PH(K,K) 
DO  211  K=1,N 

211  SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 
DO  220  1=1, N 

R(I)=R(I)+(AMU/SGL(I))*EK*GR(I) 
IF(R(I).GE.1.0)  R(I)=0. 
IF(R(I).LE.-1.0)  R(I)=0. 

220   CONTINUE 
RETURN 
END 


SUBROUTINE  PLOT(Y,N) 
DIMENSION  Y(N),X(5000) 
ISTP=N/10 
DO  10  J=1,N 
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10   X(J)=J 

c    CALL  TEK618 

c    CALL  PRTPLT(72,6) 

CALL  SHERPA ( ' PARKPARK ' , ' A ' , 3 ) 

CALL  PHYSOR(l. ,1.) 

CALL  RESET ('ALL') 

CALL  PAGE (8. 5, 11.0) 

CALL  HWROT ( ' AUTO ' ) 

CALL  XINTAX 

CALL  AREA2D(5.0,2.8) 

CALL  HEIGHT (0.10) 

CALL  COMPLX 

CALL  SHDCHR(90. 0,1, 0.002,1) 

CALL  HEAD IN ('LEARNING  CURVE$ ■ , 100 ,2 .0 , 1) 

CALL  XNAME('ITERATIONS$' ,100) 

CALL  YNAME ( ' ERR0R$ ' ,100) 

CALL  MESSAGC ADAPTIVE  LATTICE (AMU= . 5 , FIG5 . 5)$ ' , 100 , 3 . 0 , -0 . 8) 

CALL  FRAME 

CALL  GRAF(0,ISTP,N,-2.5,1.25,2.5) 
C    CALL  THKCRV(0.02) 

CALL  CURVE (X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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******Program  $********************************************************** 

*  This  program  is  designed  for  the  system  identification  experiment.    * 

*  The  LMS  algorithm  shown  in  Section  (V.C)  was  extended  to  the  linear   * 

*  phase  FIR  lattice  filter.  * 


-k 


INTEGER  ISEED,K/I,J,N,M,NB,R,MA,ML 

REAL  A(100),B(100),F(100),G(100),GD(100),YF(3000),YL(3000) 

REAL  X(100)  ,H(100)  ,YD(3000)  ,ER(3000)  ,W(3000)  ,YFF(3000)  ^LLOOOO1) 

REAL  SGL(100),PH(100,100),PS(100,100),PSD(100,100) 

REAL  GR(100),GRD(100,100),GRB(100),Y(3000)/AMU,E/C 


*  Set  Adaptation  Constant  Jl  and  Number  of  Iterations  NB 
* 

WRITE (5,1) 
1    FORMAT(5X, 'AMU' ,5X, 'NB' ) 
READ(5,*)  AMU,NB 
* 

*  Initialization 
* 

E=0. 

R=l 

M=4 

ISEED=343169 

DO  10  K=l,100 

A(K)=0 

B(K)=0 

X(K)=0 

F(K)=0 

G(K)=0 

GD(K)=0 

H(K)=0 

SGL(K)=1.0 

GR(K)=0 

GR3(K)=0 
10       CONTINUE 

DO  11  K=l.,3000 

W(K)=0 

YF(K)=0 

YL(K)=0 

YFF(K)=0 

YLL(K)=0 

YD(K)=0 

ER(K)=0 

Y(K)=0 
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11       CONTINUE 

DO  15  K=l,100 

DO  15  L=l,100 

PH(K,L)=0 

PS(K,L)=0 

PSD(K,L)=0 

GRD(K,L)=0 

15       CONTINUE 

c    C=.154 

c    A(l)=.154/C 

c    A(2)=.462/C 

C=.15 

A(1)=.15/C 

A(2)=-.45/C 

A(3)=.36/C 
* 

Random  Number  Generation 
*       (mean  zero,  unit  variance,  white  sequence) 


* 


DO   20  N=1,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 
W(N)  =  RN 

20  CONTINUE 

IF  (MOD(M,2).EQ.O)  MA=M/2 

IF  (MOD(M,2).NE.0)  MA=(M+l)/2 

IF  (MOD(M,2) .EQ.O)  ML=MA 

IF  (MOD(M/2).NE.O)  ML=MA-1 
* 

Separation  of  Coefficients 
* 

IF  (MOD(M,2).EQ.O)  THEN 
DO  21  1=1, MA 
H(I)=A(I) 
H(MA+2+I)=A(MA+l-I) 

21  CONTINUE 
H(MA+l)=A(MA+l)/2 
H(MA+2)=A(MA+l)/2 

END  IF 

IF  (M0D(M",2).NE..Q)  THEN 
DO  22  1=1, MA 
H(I)=A(I) 
H(MA+I)=A(MA+1-I) 

22  CONTINUE 
END  IF 

*  LINEAR  PHASE  FIR  FILTER 


90 


DO  30  K  =  1,NB 
X(1)=W(K) 
YF(K)=H(1)*X(1) 
IF  (MOD(M,2).EQ.O)  THEN 
DO  40  1=1, MA 

YF(K)=YF(K)+H(I+1)*X(I+1) 

40  CONTINUE 

YFF(K)=YF(K) 
DO  41  I=1,MA+1 

YFF(K)=YFF(K)+H(MA+1+I)*X(MA+I) 

41  CONTINUE 
END  IF 

IF  (MOD(M,2).NE.O)  THEN 
DO  42  1=1, MA- 1 

YF(K)=YF(K)+H(I+1)*X(I+1) 

42  CONTINUE 

YFF(K)=YF(K) 
DO  43  1=1, MA 

YFF(K)=YFF(K)+H(MA+I)*X(MA+I) 

43  CONTINUE 
END  IF 

DO  45  1=1, M 

X(M+2-I)=X(M+l-I) 
45  CONTINUE 

30    CONTINUE 

*         LATTICE  FILTER 

DO   50  K  =1,NB 

F(1)=W(K) 

G(1)=W(K) 

DO   60   I  =  1,ML 

F(I+1)=F(I)+B(I)*GD(I) 

G(I+1)=GD(I)+B(I)*F(I) 

60  CONTINUE 

YL(K)=F(ML+1) 

YD(K)=G(ML+1) 

YLL (K) =YL ( K ) +YD ( K-MA ) 

E  =  YFF(K)  -  YLL(K) 
* 

Updating  the  Reflection  Coefficients 

CALL  MLMS  (P,Q,QD,GRB,GRD,F,G,GD,B,SGL/E,AMU,ML,MA) 
ER(K)=E**2 
DO  70  J=1,ML 
GD(J)=G(J) 
70  CONTINUE 

91 


Y(K)=E 

WRITE(6,300)  K,Y(K) 

50      CONTINUE 
* 

*  Plotting  the  Learning  Curve 

CALL  PLOT(Y,N) 

300  FORMAT(3X,I6,3X,F10.5) 

STOP 

END 
* 

* 

SUBROUTINE  MLMS(PH,PS ,PSD,GRB,GRD,F,B,BD,R,SGL,EK,AMU,N,MA) 

DIMENSION  PH(100/100),PS(100,100),PSD(100,100)/GRD(100/100) 
1,GRB(100),F(100),B(100),BD(100),R(100),GR(100),SGL(101) 

GR(N)=BD(N) 

GRB(N)=F(N) 

DO  200  1=2, N 

PH(I,1)=BD(N+1-I) 

PS(I„1)=F(N+1-I) 

IT=I-1 

DO  110  K=1,IT 

PH(I,K+1)=PH(I,K)+R(N+1-I+K)*PSD(I,K) 
110   PS(I,K+1)=PSD(I,K)+R(N+1-I+K)*PH(I,K) 

DO  120  K=1,IT 
120   PSD(I,K)=PS(I,K) 
200   CONTINUE 

DO  210  K=2,N 

GR(N+1-K)=PH(K,K) 
210   GRB(N+1-K)=PS(K,K) 

DO  220  K=1,N 

DO  220  L=1,MA 
220   GRD(K,MA+2-L)=GRD(K,MA+l-L) 

DO  230  K=1,N 
230   GRD(K,1)=GRB(K) 

DO  240  K=1,N 
240   GRB(K)=GRD(K,MA+1) 

DO  250  K=1,N 
250   GR(K)=GR(K)+GR£(K) 

DO  260  K=1,N 
260   SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 

DO  270  1=1  ,N 

R(I)=R(I)+(AMU/SGL(I))*EK*GR(I) 
c    IF(R(I).GE.1.0)  R(I)=0. 
c    IF(R(I).LE.-1.0)  R(I)=0. 
270   CONTINUE 
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RETURN 
END 


* 


SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N),X(3000) 

ISTP=N/10 

DO  10  J=1,N 
10   X(J)=J 
c    CALL  TEK618 
c    CALL  PRTPLT(72,6) 

CALL  SHERPA ( ' PARKPARK ' , ' A ' , 3 ) 

CALL  PHYSOR(l. ,1.) 

CALL  RESET ('ALL1 ) 

CALL  PAGE (8. 5, 11.) 

CALL  HWROT ( ' AUTO ' ) 

CALL  XINTAX 

CALL  AREA2D(5.0,2.8) 

CALL  HEIGHT (0.10) 

CALL  COMPLX 

CALL  SHDCHR(90. 0,1, 0.002,1) 

CALL  HEAD IN ( 'LEARNING  CURVE$ '  , 100 , 2  .0 , 1 ) 

CALL  XNAME('ITERATIONS$' ,100) 

CALL  YNAME ( ' ERROR$ ' ,100) 

CALL  MESSAG( 'ADAPTIVE  G-M  LATTICE(F5. 10 ,AMU=0 . 1)$ ■ , 100 ,3 .0 , -0 . 8) 

CALL  FRAME 

CALL  GRAF(0,ISTP,N,-8.0/4.0,8.0) 
c    CALL  THKCRV(0.02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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******p  roar  am  4*************  ******  *************************************** 

*  This  program  is  designed  for  the  estimation  of  spectral  lines  in  * 

*  white  noise.  The  input  process  x(k)  consists  of  a  signal  in  noise,  * 

*  and  a  signal  may  be  a  single,  or  multiple  sinusoids.  The  algorithm  * 

*  was  derived  in  Section  (V.D).  * 
****************************************************** 

* 


* 


INTEGER  ISEED,K,I,J,N,M,NB,R,MA,ML,D 

REAL  A(100) ,B(100) ,F(100) ,G(100) ,GD(100) ,YL(5000) 

REAL  Y(100) ,YD(5000) ,W(5000) ,YLL(5000) ,YD(5000) , INP(5000) 

REAL  SGL(101),PH(100,100),PS(100,100),PSD(100,100) 

REAL  GRD(100,100),GRB(100),GR(100) 

REAL  RE(IOO) ,IM(100) ,AJ(100) ,SD,SNR,AVG,AMP, AMU,E 


*  Set  Adaptation  Constant  \l   and  Number  of  Iterations  NB 

WRITE(5,1) 

F0RMAT(5X, 'AMU' ,5X, 'NB') 

readCs,*)  AMU,NB 

* 

*  Initialization 
* 

c    ISPEC=1000 
SNR=30 . 

SD=10**(-(SNR/20)) 
AMP=SQRT(2.) 
AVG=0 . 
E=0. 
c    R=l 

PI=3. 141592654 

M=8 

MP1=M+1 

ISEED=343169 

IF  (MOD(M,2).EQ.O)  THEN 

MA=M/2 

ML=MA 

ELSE 

MA=(M+l)/2 

ML=MA-1 

END  IF 

DO  10  K=l,100 

A(K)=0 

B(K)=0 

F(K)=0 

G(K)=0 
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GD(K)=0 

SGL(K)=1.0 

GR(K)=0 

GRB(K)=0 

RE(K)=0 

IM(K)=0 

AJ(K)=0 

Y(K)=0 

10  CONTINUE 

DO  11  K=l,5000 
W(K)=0 
YL(K)=0 
YLL(K)=0 
YD(K)=0 
ER(K)=0 
INP(K)=0 

11  CONTINUE 

DO  15  K=l,100 
DO  15  L=l,100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
GRD(K,L)=0 
15       CONTINUE 

Random  Number  Generation 
(mean  zero,  unit  variance,  white  sequence) 


* 


DO   20  N=1,NB 

CALL  LNORM(ISEED,RN, 1,1,0) 
W(N)  =  SD*RN+AVG 
20   CONTINUE 

*         LATTICE  FILTER 
* 

DO   50  K  =1,NB 
AK=K-1 
c         INP(K)=AMP*COS(2*PI*.15*AK)+W(K) 
c         INP  (K)  =AMP*  ( COS  (,  2*?I *  .  1 5^AK )  +COS  ( 2*PI*  .  25*AK )  )  +W  ( K ) 

INP(K)=AMP*(COS(2^PI*.05^AK)+COS(2'VPI^.15^AK)+COS(2*?I^.25*AK 
*)+COS(2*PI*.35*AK))+W(K) 
F(1)=INP(K) 
G(1)=INP(K) 

DO   60   I  =  1,ML 

F(I+1)=F(I)+B(I)*GD(I) 
G(I+1)=GD(I)+B(I)*F(I) 
60  CONTINUE 
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* 


YL(K)=F(ML+1) 
YD(K)=G(ML+1) 
YLL(K)=YL(K)+YD(K-MA) 
E  =  YLL(K) 


*  Updating  the  Reflection  Coefficients 
* 

CALL  MLMS  (PH,PS ,PSD,GR,GRB,GRD,F,G,GD,B,SGL,E,AMU,ML,MA) 
ER(K)=E**2 
DO  70  J=1,ML 
GD(J)=G(J) 
70  CONTINUE 

C  IF  (K.NE.ISPEC)  GO  TO  50 

IF  (K.EQ.NB)  THEN 
C        WRITE(6,300)  K,INP(K),E 

c        WRITE(6,301)  K,B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9) 
c  R=R+100 

*  Determine  the  FIR  Coefficients  from  the  Lattice  Reflection 

*  Coefficients 
* 

CALL  STEPUP  (A,B,ML) 
IF  (MOD (M, 2). EQ. 0)  THEN 
DO  80  1=1 ,M/ 2 

80  A(M+2-I)=A(I) 
A(M/2+l)=2*A(M/2+l) 

ELSE 

DO  81  I=l,(M+l)/2 

81  A(M+2-I)=A(I) 

END  IF 

c    WRITE  (6,600)  (I ,A(I) , 1=1 ,MP1) 
* 

*  Calculate  the  Power  Spectrum 
* 

CALL  SPEC  (RE,IM,A,AJ,M,PI) 
D=100 
c    WRITE  (6,601)  (.005*Z,AJ(Z),Z=1,D) 


Plotting  the  Spectrum 

CALL  PLOT(AJ,D) 
C     ISPEC=ISPEC+1000 

END  IF 
50      CONTINUE 
c        WRITE  (6,300)  (YLL(K) ,K=1 ,NB) 

300  FORMAT(10(F10.7,1X)) 

301  FORMAT(1X,I5,1X,10(F10.7,1X)) 

96 


600  FORMAT  (IX, 15 , 5X,F15 .7 ) 

601  FORMAT  (IX, F5 .3 , 5X,F15 .7 ) 
STOP 

END 


SUBROUTINE  MLMS ( PH , PS , PSD , GR , GRB , GRD , F , B , BD , R , SGL , EK , AMU , N , MA ) 

DIMENSION  PH(100,100);PS(100/100),PSD(100,100),GR(100),GRB(100) 
1,GRD(100,100),F(100),B(100),BD(100),R(100),SGL(101) 

GR(N)=BD(N) 

GRB(N)=F(N) 

DO  200  1=2, N 

PH(I,1)=BD(N+1-I) 

PS(I,1)=F(N+1-I) 

IT=I-1 

DO  110  K=1,IT 

PH(I,K+1)=PH(I,K)+R(N+1-I+K)*PSD(I,K) 

PS(I,K+1)=PSD(I,K)+R(N+1-I+K)*PH(I,K) 
110   CONTINUE 

DO  120  K=1,IT 
120   ?SD(I,K)=PS(I,K) 
200   CONTINUE 

DO  210  K=2,N 

GR(N+1-K)=PH(K,K) 

GRB(N+1-K)=PS(K,K) 
210   CONTINUE 

DO  220  K=1,N 

DO  220  L=1,MA 
220   GRD(K,MA+2-L)=GRD(K,MA+l-L) 

DO  230  K=1,N 
230   GRD(X/1)=GRB(K) 

DO  240  K=1,N 
240   GRB(K)=GRD(K,MA+1) 

DO  250  K=1,N 
250   GR(K)=GR(K)+GRB(K) 

DO  260  K=1,N 
260   SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 

DO  270  1=1, N 

R( I ) =R( I ) - ( AMU/ SGL ( I ) ) *EK*GR ( I ) 

IF(R(I).GE.1.0)  R(I)=0. 

IF(R(I).LE.-1.0)  R(I)=0. 
270   CONTINUE 

RETURN 

END 
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SUBROUTINE  STEPUP  (A,B,ML) 

DIMENSION  A(100),C(100) ,B(100) 

A(l)=l. 

A(2)=B(1) 

DO  30  MINC=2,ML 

DO  10  J=1,MINC 

JB=MINC-J+1 
10   C(J)=A(JB) 

DO  20  IP=2,MINC 
20    A(IP}=A(IP)+3(MINC)*C(IP-1) 

A(MINC+1)=B(MINC) 
30    CONTINUE 

RETURN 

END 


SUBROUTINE  SPEC (RE , IM, A, AJ ,M,PI) 

REAL  RE(100),IM(100),A(100),AJ(100),Y(100) 

RE(1)=A(1) 

IM(1)=0. 

DO  91  J=l,100 

DO  92  1=1, M 

RE(I+1)=RE(I)+A(I+1)*COS(2*I*PI*.5*J/100) 

IM(I+1)=IM(I)+A(I+1)*SIN(2*I*PI*.5*J/100) 

92  CONTINUE 
AJ(J)=-10.*ALOG10(RE(M+1)*RE(M+1)+IM(M+1)*IM(M+1)) 

91    CONTINUE 

TEMP=AJ(1) 

DO  93  L=2,100 

IF  (AJ(L).GT.TEMP)  TEMP=AJ(L) 

93  CONTINUE 

DO  94  L=l,100 
AJ(L)=AJ(L)-TEMP 

94  CONTINUE 
RETURN 
END 


SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N),X(100) 
c    ISTP=N/10 
c    DO  10  J=1,N 
clO   X(J)=J 

DF=.5/N 
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X(1)=0 
DO  10  K=2,N 
10   X(K)=X(K-1)+DF 

XMIN=X(1) 
XMAX=X(N) 
XSTP=10*DF 
IYMIN=Y(1) 

IYMAX=Y(1) 

DO  20  K=2,N 

IF(Y(K).GT.IYMAX)  IYMAX=Y(K) 

IF(Y(K).LT. IYMIN)  IYMIN=Y(K) 
20   CONTINUE 

I YSTP= ( I YMAX- IYMIN ) / 5 
c    CALL  TEK618 
c    CALL  PRTPLT(72,6) 

CALL  SHERPA ( ' PARKPARK ' , ' A ■ , 3 ) 

CALL  RE SET ('ALL') 

CALL  PAGE (8. 5, 11.0) 

CALL  HWROT ( ' AUTO ' ) 

CALL  XINTAX 

CALL  AREA2D(5.0,2.3) 

CALL  HEIGHT (0.10) 

CALL  COMPLX 

CALL  SHDCHR(90. 0,1 ,0.002,1) 

CALL  HEADINC FREQUENCY  SPECTRUM$ ', 100,2 .0 , 1) 

CALL  XNAME('FREQUENCY$I ,100) 

CALL  YNAHE( 'MAGNITUDE (DB)$' ,100) 

CALL  MESSAGC FIGURE  5.21  (M=8  ,SNR=30DB)$ ' , 100 ,3 .0 , -0 .8) 
c    CALL  MESSAG( 'MODEL  ORDER  SELECTI0N(4  SINUSOIDS)$ ' , 100 ,3 . 0 , -0 .8) 

CALL  THKFRM(0.03) 

CALL  FRAME 

CALL  GRAF ( XMIN , XSTP , XMAX , IYMIN , IYSTP , IYMAX) 
C     CALL  THKCRV(0.02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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